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14.  ABSTRACT 

This  report  explores  the  mechanics  of  the  motion  (and  the  aspects  of  bio-physical  self-regulation  of  this  motion)  in  the  paramecium  s  cilium.  The  paramecium 
is  a  single-cell  animal  widely  found  in  oxygenated  aquatic  environments.  These  animals  propel  themselves,  albeit  with  limited  maneuverability,  by  the 
synchronous  motion  of  numerous  tiny  cilia  populated  around  their  flexible  bodies.  Three  theoretical  models  are  given;  (1)  a  torsional  pendulum  model  of  beat 
frequency,  (2)  a  nonlinear  seff-regulation  model  of  cilium  motion,  and  (3)  a  bio-physical  mechanism  of  hardness  control  of  the  cilium.  The  first  two  models 
reproduce  the  experimental  data  quite  accurately;  a  full  modeling  of  the  hydrolysis  of  adenosine-tri-phosphate  (ATP),  the  “currency'  of  chemical  energy  (which 
is  beyond  the  scope  of  this  report),  is  awaited  to  quantitatively  evaluate  the  third  mechanism  fully.  Essentially,  these  theoretical  results,  together  with  an 
analysis  of  the  motion  of  the  biological  cilium  and  a  comparison  of  the  results  with  various  scaled  biorobotic  hardware  models,  are  used  to  construct  a  bio¬ 
physical  description  of  a  self-regulating  mechanism  of  natural  swimming  in  the  paramecium’s  cilium  It  is  theorized  that  cross-bridge  links  between  the 
microtubule  pairs  are  the  source  of  cilium  hardness  during  the  power  stroke;  there  is  a  critical  phase  near  the  end  of  the  power  stroke  where  one  cross-bridge 
detaches  at  the  point  of  inflection  due  to  the  arrival  of  ATP,  causing  a  precipitous  reduction  in  hardness  that  signals  the  start  of  the  return  stroke,  therefore,  in 
each  beat  cycle,  there  must  be  a  reattachment  process  of  the  cross-bridge  links  and  re-hardening  of  the  cilium  during  the  early  phase  of  the  return  stroke.  Also 
discussed  is  the  invariance  of  the  elementary  oscillatory  architecture  of  integrating  sensors,  actuators,  and  controllers  from  low  to  high  Reynolds  numbers 
(paramecia  to  fish).  The  present  work  helps  lay  the  foundation  of  nonlinear  underwater  acoustic  transduction,  whose  scope  appears  to  be  uncharted. 
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EXPLORATION  OF  SELF-REGULATION  IN  THE 
NATURAL  SWIMMING  OF  THE  PARAMECIUM’S  CILIUM 


1.  INTRODUCTION 


The  paramecium  is  slipper-shaped  ciliate  protozoa  widely  found  in  oxygenated  aquatic 
environments.  Paramecia  are  100  -  350  pm  long,  deformable,  and  may  contain  up  to  3000 
flexible  cilia  populated  around  their  body,  each  cilium  being  about  1 7  pm  long  and  0.25  pm  in 
diameter  (Pemberg  and  Machemer.  1095).  For  propulsion,  each  cilium  beats  at  about  14.1  Hz  in 
water  while  undergoing  a  complicated,  three-dimensional,  phase-dependent  motion.  The 
assemblage  of  cilia  also  undergoes  a  coordinated  motion,  called  metachronic  motion — i.e.,  the 
assemblage  beats  with  a  constant  phase  difference  between  neighboring  cilia  (Machemer,  1972). 
The  paramecium  twists  while  going  forward.  It  backs  off  and  swims  in  a  new  direction  when  it 
meets  an  obstruction.  The  paramecium's  swimming  speed  is  about  two  body  lengths  per  second 
(0.44  mm/s  for  a  200-pm-long  paramecium  (Wurzel,  2003)).  The  range  of  maneuvering  motion 
that  a  paramecium  can  undergo  is  limited.  In  the  following  paragraphs,  the  literature  on  cilia 
motion  and  their  metachronic  assemblage  and  control  is  examined. 


1.1  LITERATURE  REVIEW 

Information  on  the  mathematical  treatment  of  the  three-dimensional  motion  of  a  single 
cilium  can  be  found  in  Dillon  and  Fauci  (2000),  Hines  and  Blum  (1978.  1983)  and  Brokaw 
(1966).  Further  information  on  metachronism  can  be  found  in  Blake  (1971),  Gueron  et  al. 
(1997),  Gueron  and  Levit-Guerevich  (1999a.  1999b).  Gueron  and  Liron  (1992),  Machemer 
(1972),  and  Ramia  et  al.  (1993). 

Teunis  and  Machemer  ( 1 994)  have  reported  high-speed  stereomicroscopic  measurements  of 
cilia  motion  and  have  derived  the  variation  of  curvature  and  torsion  with  the  beat  phase.  They 
have  also  identified  the  power  and  return  strokes  of  cilia  motion.  (In  the  present  work,  the 
Teunis  and  Machemer  results  are  revisited.) 

Dute  and  Kung  (1978)  have  given  the  physiology  of  a  cilium,  indicating  how  the  cilium 
motion  is  mechanically  produced  when  a  torque  is  applied  at  the  cilium's  base.  Cilium  motion  is 
generated  by  a  combination  of  a  pair  of  linear  tubular  actuators  (microtubules)  at  the  center, 
which  are  surrounded  by  nine  pairs  of  push-pull  microtubules  spaced  at  40°  apart  (this  actuator 
configuration  of  microtubules  is  called  the  “(9  +  2)  structure").  Dute  and  Kung  offer  a 
mechanical  model  that  bears  some  similarity  to  helicopter  swashplates. 

Hines  and  Blum  (1978.  1983)  have  carried  out  a  nonlinear  theoretical  analysis  of  the  motion 
of  a  flagellum  (approximately  a  longer  cilium)  using  a  sliding  filament  model.  They  reproduce 
bending  propagation  and  argue  that  the  cilium  must  be  twisting.  Dillon  and  Fauci  (2000)  have 
hydrodynamically  modeled  the  cilium,  incorporating  a  most  detailed  (9+  2  “axenome")  internal 
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elastic  actuation  structure.  In  their  study,  they  reproduce  the  cilium  beating  as  an  emergent 
property  of  the  coupling  of  the  fluid  and  the  actuator  architecture.  Gueron  and  Levit-Guerevich 
(2001 )  have  incorporated  the  9  +  2  architecture  in  their  theoretical  model. 

Cilium-cilium  hydrodynamic  interaction  has  been  theoretically  examined  by  Gueron  and 
Liron  (1992)  and  Gueron  et  al.  (1997),  extending  the  substantial  earlier  theoretical  work  on  the 
hydrodynamics  of  a  single  cilium  (see  their  paper  for  the  earlier  works).  Their  model  is  two- 
dimensional.  They  show  that  in  two  cilia  beating  at  a  random  phase  difference,  they  synchronize 
within  two  beat  cycles  when  brought  closer  together;  that  is,  the  metachronal  pattern  develops 
autonomously.  The  implication  is  that  metachronism  is  a  self-organizing  phenomenon. 

How  the  motion  of  the  cilium  (that  is,  the  hydrodynamics)  is  controlled  can  be  gleaned  from 
the  following.  Machemer  (1972),  Nakaoka  et  al.  (1984),  and  dePeyer  and  Machemer  (1983) 
have  shown  that  voltage  stimulation  affects  cilia  motility.  Slow  voltage  ramps  (5  to  10  mV/s) 
and  moderate  amplitudes  (±20  mV)  are  shown  to  make  the  membrane  properties  time  dependent, 
which  affects  ciliary'  response  (cycle  time  and  relaxation  time).  Pemberg  and  Machemer  (1995) 
and  Deitmer  et  al.  (1984)  have  shown  that  the  cilia  beating  frequency  and  reverse/forward 
motion  switching  are  affected  by  calcium  channel  activation.  Jaffe  (2007)  has  examined  60  data 
sets  on  flagella  and  cilia  and  has  shown  that  stretch-activated  calcium  channels  cause  calcium 
and  other  waves  to  propagate  at  speeds  of  100  to  1000  pm/s.  which  affects  the  tubular  actuators 
mentioned  earlier. 

Understanding  of  the  cilium  bending  at  the  cellular  level  due  to  cross-bridge  linking  has 
made  rapid  strides.  Cordova  et  al  (1992)  have  mechanically  modeled  the  attachment  of  a  single 
motor  molecule.  Electrostatic  force  between  the  cross-bridge  and  the  binding  site  is  used  to 
generate  the  cross-bridge  link  motion,  and  against  an  elastic  restoring  force,  the  cross-bridge 
fluctuates  and  its  step  size  varies  nonlinearly  with  force. 

Tamm  and  Tamm  (1989)  have  studied  microtubule  sliding  experimentally;  the  attachment  of 
a  Ca  ion  to  a  cilium  is  found  to  involve  locally  confined  oscillatory  bending,. 

Some  say  that  the  mechanism  of  metachronism  is  hydrodynamic  coupling,  and  some  say 
that  it  is  electrical  (membrane  potential)  in  origin  (Tamm,  1984).  The  former  is  largely  based  on 
two-dimensional  fluid  dynamics  modeling,  and  the  paramecium's  cilia  also  show  electro¬ 
mechanical  response.  (The  present  report  explores  these  approaches  from  the  point  of  view  of 
self-regulation  in  mechanical  systems.) 

A  photon  correlation  technique  has  led  to  the  measurement  of  relaxation  times  of  5  ps  in 
single  muscle  fiber  (Yeh  et  al.,  1990).  A  lower  frequency  of  370  kHz  is  also  present.  The  high 
frequency  comes  from  the  cross-bridge  linking  and  the  low  frequency  from  certain  thick 
filaments. 

Bouzarth  et  al.  (2007)  have  carried  out  modeling  and  experiments  on  the  submicron-scale 
flow  of  cilia.  They  found  a  strong  coupling  between  cilium  motion  and  the  flow:  the  fluid  orbit 
is  epicyclical  (related  to  the  slow  motion  of  the  cilium)  with  coherent  fluctuations  (related  to  the 
precession  rate). 
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Guirao  et  al.  (2010)  have  shown  that  the  alignment  of  ciliary  beating  to  fluid  flow  direction 
is  controlled  by  a  mechanism  in  which  certain  cilia  act  as  sensors  for  lowpass  filtering  of 
hydrodynamic  forces  and  for  generating  a  polarity  signal  for  directional  alignment. 

Small  biorobotic  cilia  have  also  been  fabricated,  although  their  hardness  remains  unchanged. 
Kongthon  et  al.  (2010)  have  shown  that  cilia-based  PDMS  (polydimethylsiloxane)  actuators  have 
lower  resonant  frequency  in  liquid  than  in  air.  They  show  that  this  can  be  explained  by  added 
mass  effects  only.  Shields  et  al.  (2010)  describe  an  experiment  on  PDMS  biomimetic  cilia 
stirring.  The  Peclet  number,  which  is  a  ratio  of  advective  to  diffusive  motion,  is  estimated  to  be 
1 0;  that  is,  their  cilium  flow  is  advection  dominated.  In  both  of  these  papers,  phase  dependent 
bending  is  absent.  They  find  the  flow  within  the  conical  cilia  orbit  to  be  non-unidirectional,  but 
unidirectional  beyond  the  tip,  and  they  term  it  “mixing'’  and  “pumping,”  respectively.  They 
model  the  flow  to  be  a  linear  superposition  of  shear  and  pressure-driven  flow  between  parallel 
plates. 

Theoretical  formulation  of  the  equations  of  motion  of  a  cilium  (and  flagellum)  remains  a 
difficult  problem.  Gueron  and  Liron  (1992)  have  applied  corrections  to  Lighthill's  theory. 
Although  a  large  body  of  analytical  research  on  flagellum  motion  was  carried  out  during  the 
1970s,  that  research  (and  interpretation  of  the  research's  measurements)  did  not  predict  that  the 
motion  is.  in  fact,  three-dimensional.  The  finding  by  Hines  and  Blum  (1978.  1983)  that  there  is 
twist  indicates  an  improvement  in  our  understanding  of  how  three-dimensionality  evolves  in  the 
cilium  actuator.  Robotic  modeling  may  be  viewed  as  a  supplementary'  tool  to  our  understanding 
of  three-dimensionality.  Therefore,  it  would  be  useful  to  build  a  robotic  model  of  the  cilium 
where  the  phase-dependent  deformations  are  produced  by  summing  independent  orthogonal 
oscillatory  motions  (lOOMs),  as  in  typical  coordinate  systems  (robotic  simulations  using  coupled 
actuators  make  programming  difficult).  (The  present  report  considers  what  these  IOOMs  are  and 
explores  the  origin  of  coupled  motions.  Contemporary  PDMS  cilia  actuators  do  not  accurately 
reproduce  the  cilium  motion,  and  it  is  hoped  that  the  IOOMs  will  provide  a  solution.) 

The  cilium  topic  is  of  interest  because  it  is  a  model  of  microscopic,  elementary,  natural 
swimming  propulsion.  A  substantial  amount  of  experimental  and  theoretical  information  is 
available  on  its  physiology,  motion,  theoretical  hydrodynamic  modeling,  metachronism,  and 
control.  Metachronism  by  itself  is  a  very  interesting  example  of  cilium-cilium  coupling.  There 
is  evidence  of  self-organization,  as  well  as  evidence  of  ionic  control  of  cilium  beating.  The 
potential  relationship  of  natural  swimming  and  self-organization  plus  ionic  control  has  far- 
reaching  consequences  in  the  sense  that  hydrodynamics  and  control  may  have  co-evolved — see 
Llinas  (2001). 


1.2  PURPOSE  OF  THIS  WORK 

The  present  work  addresses  autonomy  and  three-dimensionality.  The  focus  is  on  different 
sources  of  resonant  oscillations.  Also  considered  is  how  three-dimensionality  develops  in  a  solid 
surface,  and  how  the  development  of  this  three-dimensionality  and  control  are  integrated.  The 
approach  taken  is  that  modeling  of  how  the  different  aspects  of  the  problem  are  integrated  as  a 
whole  is  a  fruitful  way  of  understanding  autonomy. 
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The  degree  of  three-dimensionality  and  control  in  the  cilium  varies  with  the  power  and 
return  strokes.  Hence,  it  is  important  to  determine  their  phase  accurately.  How  long  is  the 
propulsion  stroke?  Per  Teunis  and  Machemer  (1994),  it  is  1/7  to  2/7  of  the  beat  period.  This 
finding  is  revisited  here  because  the  cilium  motion  shows  that  near  the  extremities  of  the  power 
and  return  strokes,  the  outer  and  inner  parts  of  the  cilium  change  the  stroke  at  different  phases, 
which  makes  it  difficult  to  determine  the  exact  ends  of  the  strokes.  For  self-regulation  to  be 
explored,  the  phase  transition  between  the  power  and  return  strokes  needs  to  be  determined 
accurately. 

Voltage  stimulation  experiments  have  shown  that  ion  flow  in  the  microtubules  controls  the 
torque  produced.  Theoretical  work  has  shown  that  (1)  the  cilium  motion  emerges  when 
appropriate  moments  are  applied  at  the  cilium  microtubules,  and  (2)  neighboring  cilia  interact 
and  phased  motion  ensues  between  them  autonomously.  The  links  between  autonomous 
regulation,  ion  flow-based  control,  and  cilium  motion  need  to  be  examined  because  they  may 
provide  clues  as  to  how  hydrodynamics  and  control  are  integrated  in  natural  swimming. 

To  understand  the  motion  of  a  cilium.  this  report  brings  in  several  tools  that  have  not  been 
used  previously  in  cilium  investigations.  These  tools — Frenet’s  equation,  nonlinear  oscillatory 
modeling,  biorobotics,  and  materials  science — provide  an  opportunity  to  integrate  structure, 
hydrodynamics,  and  control,  whereby  a  more  internally  consistent  understanding  can  be  reached, 
and  perhaps  applied  to  nonlinear  transduction,  which  is  not  well  developed.  Our  physical 
rendition  is  a  way  to  understand  internal  consistency  when  different  disciplines  such  as 
hydrodynamics  and  control  are  being  integrated. 

In  the  following  sections,  a  theoretical  model  on  torsional  oscillation  of  the  cilium  is  first 
developed,  and  this  is  followed  by  a  nonlinear  model  of  the  moments  applied  at  the  base.  After 
this,  the  motion  of  the  biological  cilium  is  examined  by  revisiting  Teunis  and  Machemer' s 
(1994)  cilium  orbits.  The  relationship  of  the  analysis  of  the  cilium  motion  to  the  models  is 
examined.  Cilium  motion  is  robotically  reproduced  in  a  scaled-up  model  and  compared  with  that 
of  the  biological  cilium.  Finally,  a  synthesis  of  the  results  is  given  and  a  modified  fracture 
theory  is  applied  to  explain  the  origin  of  variation  in  hardness.  The  process  of  self-regulation  is 
discussed.  The  report  contains  a  substantial  amount  of  supporting  information  in  the  appendix, 
which  covers  the  following  topics: 

•  Two-dimensional  hydrodynamic  modeling  based  on  Frenet's  equations  assuming  the 

cilium  to  be  composed  of  hinged  links. 

•  Robotic  exploration  of  the  independent  orthogonal  oscillatory  motions  that  produce  the 

cilium  motion. 

•  Development  of  “frictionless"  actuators,  continuing  the  theme  of  resonant  oscillation 

when  one  expects  the  amplitude  of  oscillation  and  quality  factor  to  increase. 

•  Hardware  operation  of  the  cilium  IOOMs  via  analog  olivo-cerebellar  dynamics. 

•  Experimental  investigation  to  understand  how  torsional  modulus  affects  the  development 

of  three-dimensionality. 

•  Consideration  of  how  the  microcosm  of  actuation  and  control  in  the  tiny  paramecium's 

cilium  could  conceivably  be  related  to  natural  swimmers  at  higher  Reynolds  numbers. 
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2.  MODELING 


Here,  natural  swimming  is  defined  as  one  where  the  cilium  actuator  is  forced  at  the  natural 
frequency  of  oscillation  and  the  resulting  motion  is  describable  by  dynamic  systems,  is  nonlinear 
with  the  potential  for  chaotic  motion,  and  has  a  self-regulating  nature. 


2.1  MODELING  OF  THE  NATURAL  FREQUENCY  OF  OSCILLATION  OF  THE 
BIOLOGICAL  CILIUM 

Consider  the  cilium  to  be  a  torsional  pendulum  (as  shown  in  figure  la)  subjected  to 
opposing  torque  r  and  drag  D.  For  small  rotational  angles  0 , 

T  =  -kO, 

where  T  k  is  a  torsional  spring  constant.  If  /  is  the  moment  of  inertia, 

w  =  r  =  -Tke, 

and  the  solution  is,  for  ru  =  2 with  fa  being  the  natural  beat  frequency  (1/s), 
ru  —  yJT k/  /  . 

In  our  case,  for  Reynolds  number  Re  ~  0.1  and  drag  coefficient  Cp,  CD  «  Re  (/  is  constant), 

ru  oc  yp~k  OC  yjCn  °C  y[Re  , 
which  gives 
/,  «  Vto. 

For  the  same  cilium  in  different  fluids. 


where  u  is  the  kinematic  viscosity,  /r  is  the  dynamic  viscosity,  and  p  is  the  density  of  the  fluid. 
Assuming  p  to  be  constant, 

f„=K( Vyfc), 
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where  K  is  a  constant.  Therefore,  as  p  drops,  the  oscillation  frequency  f0  increases.  This 
relationship  is  compared  with  Machemer’s  measurements  in  figure  lb.  The  trend  is  in  agreement 
with  Machemef  s  (1972)  measurements.  In  view  of  the  variable  hardness  of  the  cilium  during 
the  power  and  return  strokes,  this  relationship  is  approximate. 


(a)  (b) 

Figure  1.  (a)  Schematic  of  the  Torsional  Pendulum  Model 
of  the  Cilium;  (b)  Comparison  of  the  Model 
(solid  line:  f0  =  K(l/Vp),  K  =  75,  and  p  in  cp)  with 
Measurements  (dots)  Due  to  Machemer  (1972)) 


2.2  DYNAMIC  SYSTEMS  MODELING  OF  CILIUM  MOTION 

In  this  section,  olivo-cerebellar  dynamics,  which  owes  its  origin  to  the  inferior-olive  (10) 
neuron  dynamics  work  of  Llinas  and  Yarom  (1981),  and  Kazantsev  et  al.  (2003,  2004),  is  applied 
to  the  cilium  motion.  Direct  evidence  of  the  existence  of  such  neurons  in  paramecium  is  lacking. 
However,  there  is  corroborating  evidence  of  a  role  of  Ca  ions  in  paramecium  and  cilium  motion 
control  (describable  probably  as  a  primitive  neuron),  and  one  can  assume  that  each  cilium  has  its 
own  ion-related  control.  At  the  very  least,  it  is  noted  that  IO  neurons  are  merely  one  of  many 
situations  in  nature  where  equations  such  as  that  in  equation  ( 1 )  apply. 

Following  Kazantsev  et  al.  (2003),  the  model  of  the  /h  ion-related  controller  is  given  as 
follows: 
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(1) 


z, 

P,:(Z,)-W, 

'0 

— 

+ 

f(  a(Zi  ~  ha)_ 

_~eCa. 

where  the  variables  z,  and  wt  are  associated  with  sub-threshold  oscillations  and  low-threshold 
(G/-dependent)  spiking  in  olivo-cerebellar  dynamics  (the  higher  threshold  (Na  -dependent) 
spiking  oscillator  of  Kazantsev  et  al.'s  model  is  not  considered).  The  constant  parameter  Eca 
controls  the  oscillation  time  scale:  and  Ica  drives  the  depolarization  levels.  The  nonlinear 
function  is 


pl:(z,)  =  zl(zl  -  a, )(1  -z(). 

The  function  Iexl  ( t )  is  the  extracellular  stimulus,  whose  amplitude  and  duration  are  used  for 
the  purpose  of  control  (changing  the  motion  of  paramecium  direction,  for  example).  If  Iext  (/)  =  0 
(independent  oscillator),  the  nonlinear  function  is  given  by  pr  ( z, )  =  z,  (z,  -  a ,  )(1  -  z, ) ,  where  a , 
is  a  constant  parameter  associated  with  the  nonlinear  function.  Equation  (1)  can  be  written  as 

z,  +F(z,)zl  +  kz,  +sl  =  0,  (2) 

where  F  is  a  cubic  polynomial  function  and  k  is  a  constant.  Equation  (2)  bears  resemblance  to 
Lienard's  oscillator  (in  contrast,  the  function  F  is  a  well-defined  quadratic  in  the  familiar  van  der 
Pol's  oscillator)  (Khalil,  1996;  Slotine  and  Li,  1991).  The  oscillator  exhibits  a  closed-orbit  E,  in 

the  state  space  (  z,  —  z, );  that  is,  (  z,  —  wt ),  which  is  also  known  as  limit  cycle  oscillation  (LCO), 
the  constant  parameters  determining  the  form  of  E, .  Equation  (1)  is  solved  using  the  analog 
oscillators  described  in  Bandyopadhyay  et  al.  (2008). 

The  (z,  w)  waveform  data  are  gathered  from  Simulink  and  saved  into  files  for  reading;  these 
are  z,  w  waveforms,  their  derivatives,  and  their  second  derivatives  (not  shown  in  figure  2  due  to 
the  presence  of  an  intlection  point  in  the  first  derivative).  In  figure  2,  the  phase  maps  are  scaled 
to  the  same  size  as  the  experimental  data  and  centered  at  the  same  point  so  that  waveforms  can 
be  compared.  The  modeled  phase  maps  are  also  rotated  to  the  proper  orientation,  which  is  a  10° 
rotation  for  position  data  (z,  w)  and  -55°  for  derivative  data  (z,  vv). 

The  modeled  states  (z,  w)  are  compared  with  measurements  in  figure  2.  Note  that  the 
experimental  data  are  for  axes  shown  in  figure  13  for  the  biological  cilium  for  a  given  value  of  s. 
The  agreement,  even  in  the  first  derivative,  is  rather  good  for  a,  =  0.01 5.  The  simulation  shows 
that  the  cilium  follows  a  track  that  is  descrihahle  as  a  limit  cycle.  Therefore,  it  has  an  autonomic 
character ;  the  cilium  is  controlled  without  sensors. 

The  (i  =  0,  vv  =  0)  location  is  an  inflection  point  (figure  2b)  in  both  the  biological  cilium 
and  the  modeled  state  map.  This  is  a  new'  finding  and  that  is  significant  from  the  point  of  view 
of  the  nonlinear  oscillatory'  behavior  of  the  cilium. 
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(a)  (b) 

Figure  2.  Comparison  of  the  Tracks  of  the  Cilium  of  a  Paramecium  as  Measured  by 
Machemer  (1972)  and  Per  the  Present  Authors'  Computation  of  LCO  (Limit  Cycle 
Oscillation)  Using  the  Lienard  Oscillator  (Equation  (2)) 

(The  model  is  in  blue,  and  the  cilium  distal  point  data  are  in  red.  Position  is  plotted  on  the  left 
and  velocity  on  the  right.  The  value  of  a,  is  0.015.  (a)  (z,  w),  and  (b)  z.w.  Arbitrary  units. 


The  ( z ,  vi’)  waveforms  and  the  attractor  basins  have  been  computed  (not  shown  here).  There 
is  an  option  to  show  the  component  phases  of  each  waveform.  The  user  may  also  show  a  certain 
percentage  of  the  data  saved  by  Simulink,  with  small  values  showing  the  limit  cycle  and  large 
values  showing  the  entire  limit  cycle  approach  from  the  initial  condition.  For  comparing 
waveforms,  a  small  percentage  (20%  or  lower)  is  most  convenient.  To  see  the  entire  attractor 
basin  paths  from  multiple  initial  conditions,  the  user  can  show  100%  of  the  data,  which  is  a  time- 
consuming  process.  These  graphs  may  be  obtained  from  the  lead  author.  When  the  paramecium 
hits  an  obstacle,  its  return  to  natural  beating  would  be  described  by  these  attractor  maps 
(4,*0). 
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3.  BIOLOGICAL  AND  ROBOTIC  CILIUM  DATA 


A  rigid  link  modeling  of  the  hydrodynamics  of  the  cilium  has  been  undertaken  by  our  co¬ 
worker  Norman  Toplosky  (NUWC  Code  1 532)  (see  section  A.l  in  the  appendix).  Thus  far,  a 
two-dimensional  simulation  has  been  completed,  and  the  three-dimensional  simulation  results 
are  awaited.  The  cilium  stiffness  is  arbitrarily  varied  between  the  power  and  the  return  strokes. 
An  asymmetric  beating  pattern  is  produced  (figure  A-2).  The  base  force  integral  over  the  beat 
cycle  shows  some  positive  thrust  (figure  A-3).  In  the  three-dimensional  simulation,  higher  thrust 
is  expected  to  be  produced.  To  complement  this  work,  the  present  authors  undertook  the 
approach  of  biorobotic  simulation  of  the  fluid-loaded  cilium  and  an  analysis  of  the  biological 
cilium.  drawing  from  the  published  data  of  Teunis  and  Machemer  (1994). 


3.1  REVISITING  BIOLOGICAL  CILIUM  DATA 

3.1.1  Collection  and  Processing  of  the  Tennis  and  Machemer  (1994)  Data 

To  revisit  Teunis  and  Machemer's  (1994)  experiments,  their  orthogonal  cilium  stereo¬ 
microscopy  pictures  were  electronically  digitized  in  pixels;  the  cilium  motion  given  by  the 
position  vector  (in  pixels)  of  a  space  curve  r(sj)  =  {x(.v,/),  y(sj),  z(s,t)}.  This  is  shown  for  the 
propulsion  longitudinal  plane  in  figure  A-4  of  the  appendix.  The  cilium  is  divided  into  20 
segments  (initially  of  non-uniform  spacing).  Define  s  as  the  axial  length  from  the  base  to  the 
distal  point  ( L )  of  the  particular  cilium  segment;  also  define  an  orthogonal  coordinate  system  (x. 
y ,  z)  that  represents  the  propulsive  (or  stream-wise),  vertical,  and  spanwise  directions, 
respectively,  and  t  is  time. 

The  University  of  British  Columbia  (UBC)  Biology  Department  has  produced  a  useful  gif 
animation  of  the  cilium  motion  using  these  stereoscopic  data  (see  figure  13a.  left  column,  326- 
kB  gr/image).  In  this  report,  the  term  “revisited  TM  data"  pertains  to  our  (pixel)  digitization  of 
the  cilium  motion  as  it  appears  in  the  UBC  gif  graph,  with  the  data  due  to  Teunis  and  Machemer 
(1974),  and  their  subsequent  parametric  processing  and  smoothing.  The  graphs  used  for 
digitization  are  from  the  UBC  gif  file  because  the  cilium  position,  velocity,  and  acceleration  do 
not  visually  exhibit  any  spatio-temporal  kink  (smooth  gradients  of  time  and  length). 

Due  to  cilium  bending,  initially,  the  positions  (in  pixels)  along  the  length  were  not  equally 
spaced  in  the  longitudinal  and  top-view  planes.  For  this  reason,  the  digitization  had  to  be 
processed  in  the  following  manner.  With  cilium  bending,  the  (y(s,t))  distances  could  be 
determined  more  accurately  from  the  (x(t).y(t))  view  (side  view),  and  the  ( x(s,t ),  z(s,t))  position 
values  could  be  determined  from  the  (x(t),  z(t))  view  (top  view).  The  initial  length  was  obtained 
from  the  three  coordinates  obtained  this  way.  From  the  polynomial  fits  at  each  time  step  (TS), 
the  cilium  was  digitized  at  equal  length  intervals  (r(s,  /)).  The  data  were  further  smoothed  along 
the  orbits.  Subsequently,  a  three-dimensional  parametric  fit  (figure  3c)  was  applied;  see 
Zwillinger  (1996).  The  result  was  a  smoothed  position  distribution  both  along  the  cilium  lengths 
at  each  time  step  and  also  of  their  rotational  orbits  obtained  from  position  vectors  at  equispaced  s 
locations. 
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The  improvement  in  the  accuracy  due  to  axial  and  rotational  track  smoothing  of  the  cilium 
position  vector  data  is  given  in  the  tables  A-2  through  A-4  in  the  appendix.  The  cilium  length  of 
17  pm  is  maintained  at  all  time  steps  (the  diameter  is  0.25  pm)  (Pemberg  and  Machemer,  1995); 
the  local  lengths  are  also  maintained  at  all  positions  along  the  length  of  the  cilium  at  all  time 
steps.  Figure  6  of  Machemer  (1972)  shows  that  the  biological  cilium  data  used  here  are  for  water 
(viscosity  of  40  cP).  Figure  3  show  s  the  variation  of  the  cilium  surface  length  for  one  beat  cycle, 
dividing  the  cilium  into  10  lengths;  the  top  plot  shows  the  raw'  data,  the  middle  plot  shows  the 
smoothed  data  and  the  bottom  plot  shows  the  parametric  fit  to  the  data  smoothed  along  length 
and  orbit. 

Top  and  side  views  of  biological  cilium  data,  with  time  step  (TS)  numbers,  are  shown  in 
figure  4.  These  plots  show7  that  the  beat  period  is  composed  of  two  kinds  of  cilium  forms;  the 
cilium  remains  fairly  straight  from  TS  38  to  14  and  its  form  is  substantially  curved  from  TS  14  to 
38.  The  cilium  would  impart  greater  momentum  into  the  fluid  when  it  produces  higher  velocities 
in  the  outer  region;  therefore.  TS  38  to  14  are  likely  to  indicate  the  power  stroke.  During  the 
power  stroke,  the  cilium  stays  fairly  close  to  the  (z  =  0)  plane.  Therefore,  the  cilium  trajectory  is 
shaped  approximately  like  the  capital  letter  D.  This  observation  is  used  for  the  robotic 
simulation;  see  section  A. 3  in  the  appendix. 
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3D  Length(s)  of  Chosen  Path(s) 
Biological  Cilia -17  um  length  -07 -Dec-2011 
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Parametric  Fit  3D  Length(s)  of  Chosen  Path(s) 
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Figure  3.  Variation  of  the  Ciliunt  Axial  Distance  from  the  Base  with  Beat  Phase  (=  Time  Step 
Number)  in  the  Digitized  Biological  Cilium  in  the  Revisited  TM  Data — (a):  Raw  Distribution; 
(b):  After  Smoothing  Along  Cilium  Length;  (c):  After  Parametric  Fit  of  the  Closed  Orbits  to 

Smoothed  Axial  Data 
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Figure  4a.  Top  View  (x,  z)  of  the  C ilium  Beat  Cycle  in  the  Revisited  TM  Data  (Axes:  pm.) 

3D  plot  of  fitted  data 

Bio  Cilia  Data  -  17  um  length  -  10-Sep-2010 
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3D  plot  of  fitted  data 

Bio  Cilia  Data  -  17  um  length  -  10-Sep-2010 
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Figure  4b.  Side  View  (x,  y)  of  the  Cilium  in  the  Plane  of  the  Power  Stroke  During  the 
Beat  Cycle  in  the  Revisited  TM  Data  (Axes:  pm.) 
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3.1.2  Analysis  of  the  Teunis  and  Machenter  (1994)  Data 

Using  the  smoothed  biological  cilium  data,  the  velocities  and  accelerations  in  the  ( x{t).y{t ), 
z(t))  directions  were  calculated.  The  units  are  as  follows:  length  (pm),  velocity  (pm/s), 
acceleration  (gm/s2),  curvature  (rad/gm),  and  torsion  (rad/gm). 

The  straightness  of  the  cilium.  which  is  a  measure  of  its  hardness,  is  considered  in  figure  5. 
The  cilium  is  divided  into  20  axial  segments.  Figure  5  shows  the  distributions  of  (a)  the  angle  of 
the  velocity  from  the  x-axis,  and  (b)  the  angle  between  the  axial  length  segments  and  the  x-axis. 
See  section  A.2.5  for  definitions.  In  (a),  the  velocity  angles  are  nearly  equal  in  all  length 
segments  from  45°  to  0°  to  90°  values,  a  duration  that  is  identified  as  the  power  stroke.  On  the 
other  hand,  during  the  return  stroke,  the  velocity  angles  vary  widely  from  90°  to  180°  to  45°. 

The  velocity  angles  along  the  cilium  length  fan  out.  pivoting  near  TS  13  in  (a).  The  significance 
of  this  phase  is  considered  later.  In  (b),  the  cilium  position  angle  also  varies  w  idely  during  the 
return  stroke. 

The  velocities  are  shown  in  figure  6,  and  the  accelerations  are  shown  in  figures  7  and  8. 
During  the  power  and  return  strokes,  the  velocities  lie  approximately  in  two  separate  orthogonal 
planes.  The  velocities  are  high  during  the  power  stroke  and  are  reduced  during  the  return  stroke. 
Shields  et  al.  (2010)  report  tip  speeds  in  their  rigid  biological  cilium  of  700-800  pm/s  and  tracer 
speeds  of  200  pm/s.  Figure  7  shows  that  the  cilium  produces  an  inward  ’’tornado.”  with  the 
magnitude  increasing  away  from  the  base.  Two-dimensional  cilium  motion  cannot  produce  the 
tornado.  The  acceleration  limits  drop  by  a  factor  of  3  if  the  cilium  length  is  reduced  from  50  to 
17  pm. 

The  velocity  field  folds  at  (i  =  0,y  =  0,z  =  0).  The  time  instant  occurs  nearly  at  the  shift 
from  the  power  stroke  to  the  return  stroke.  This  time  instant  appears  in  figure  2b  as  a  point  of 
inflection.  As  per  the  nonlinear  modeling  given  earlier,  this  phase  denotes  neutral  equilibrium  in 
a  nonlinear  pendulum  (this  is  the  vertical  position  in  a  nonlinear  pendulum  if  the  string  is 
replaced  by  a  solid  rod;  the  topology  is  a  saddle — a  neutral  equilibrium  position  that  is  very 
sensitive  to  disturbances;  1 80"  later,  it  is  a  focus — a  stable  phase).  The  physical  significance  of 
this  phase  is  further  explored  later. 
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Cilium  velocity  angle  from  x-axis  for  Chosen  Path(s) 
Biological  Cilia  •  17  um  length  -  10-Jan-2012 


Step  number 


Cilium  angle  from  x-axis  for  Chosen  Path(s) 
Biological  Cilia  - 17  um  length  -  10-Jan-2012 


Figure  5.  Variation  of  the  Local  Axial  Angle  of  the  Orbital  Velocity  and  Position 
of  the  Biological  Cilium  Length  Segments  with  Respect  to  the  x-Axis  (Power 
Stroke  Axis)  During  the  Beat  Cycle  in  the  Revisited  TM  Data 
(The  graph  is  independent  of  cilium  length  (50  or  17  pm).  Color  code:  Different 
color  signifies  every  5%  of  cilium  length  starting  from  the  base.  The  horizontal 
lines  in  (a)  are  at  45°  and  90°  angles.) 
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Smoothed  Xdot  vs.  Ydot  vs.  Zdot 
Biological  Cilia  - 17  um  length  -07 -Dec-2011 
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Figure  6.  Variation  ofXdot(t)  (=  x(t))  vs  ydot(t)  (=  y(t))  vs  zdot(t)  (=  z(t)) 

Using  the  Revisited  TM  Data 
(Color  code:  Different  color  every  10%  of  cilium  length 
starting  from  the  base;  velocities  are  in  fim/s.) 


Smoothed  X,  Y,  Z  Accelerations 
Biological  Cilia  - 17  um  length  •  07-Dec-2011 


3000 


0 

"5T 

I  * 

IQ 

«  -10 


V  (vertical) 


x  (propulsion) 


Figure  7.  Inward  Acceleration  for  the  Revisited  TM  Data 
(Color  code:  Red  loops  of  constant  cilium  length  appear  every  10%  of 
cilium  length;  acceleration  arrows  are  in  blue;  acceleration  is  in  pm/s2.) 


Average  X,  Y,  Z  Accelerations  Along  Length 


Phase 

Figure  8.  Calculated  Acceleration  Components  (pm/s2)  in  the 
Biological  CiUum  Using  the  Revisited  TM  Data 


How  long  is  the  power  stroke?  There  is  some  ambiguity  as  to  the  boundary  of  the  power 
stroke.  Teunis  and  Machemer  (1994)  do  not  explicitly  define  the  extent  of  the  power  stroke  but 
seem  to  suggest  it  to  be  from  where  .v-acceleration  crosses  the  zero  value  to  where  it  reaches  the 
maximum  value.  Defining  the  power  stroke  as  when  the  cilium  is  moving  in  the  negative 
propulsion  direction  seems  most  reasonable,  thereby  pushing  the  cilium  in  the  positive 
propulsion  direction.  This  is  also  when  acceleration  is  increasing  in  figure  8 — approximately 
from  TS  35  to  12.  (The  end  of  the  power  stroke  can  be  ambiguous  because  the  proximal  part 
starts  the  return  stroke  prior  to  the  distal  part.).  The  total  duration  of  the  power  stroke  is  0.50  of 
the  beat  period,  which  is  larger  than  the  range  (0.14  to  0.28)  arrived  at  by  Teunis  and  Machemer 
(1994).  The  extent  of  the  power  stroke  cannot  be  accurately  determined  from  the  cilium  position 
data  alone;  accurate  calculation  of  acceleration  is  a  more  reliable  indicator. 

An  average  radius  of  the  cilium  Rav(,  was  defined,  as  given  in  section  A.2  of  the  appendix. 
The  motion  of  the  cilium  is  modeled  as  a  disc  whose  velocity  varies  radially  (Wakeling  and 
Ellington,  1997).  The  disc  area  is  divided  equally  at  Rm,g  such  that  the  momentum  imparted  to 
the  fluid  by  the  inner  and  the  outer  areas  is  equal.  The  variation  of  R^g  with  phase  is  shown  in 
figure  9. 
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Ravg  in  y,z  plane  of  Chosen  Path(s) 
Biological  Cilia  - 17  urn  length  -07-Dec-2011 


Figure  9.  Variation  of  Calculated  Average  Radius  (Ravg  in  uni)  of  the  Biological  Cilium 
During  the  Beat  Cycle  in  the  Revisited  TM  Data 
(Location  P  is  a  point  of  inflection  on  the  cilium  at  TS  14;  see  figure  1 7a.) 


The  RUVg  (this  is  a  curved  length  radius  average)  values  were  multiplied  by  the  local  velocity 
in  the  propulsion  direction  (x-axis)  to  determine  the  directionality  and  relative  scale  of  the 
momentum  imparted  to  the  fluid.  The  results  are  shown  in  figure  10.  At  TS  15.  R<ng  is  a 
minimum:  P  is  a  point  of  inflection  on  the  cilium  at  TS  14.  just  prior  to  when  Ravg  is  a  minimum. 
The  variable  ( velocity  x  Rayg  length)  crosses  zero  at  TS  35  and  14  (figure  10).  This  is  an 
indication  that  the  power  stroke  extends  from  TS  38  to  14.  and  the  return  stroke  extends  from  TS 
14  to  38.  Recall  that  the  demarcating  phases  (TS  14  and  38)  have  been  identified  as  saddle  and 
focus,  respectively.  Figures  9  and  10  show  that  the  reduction  of  R„yg  (the  reduction  in  the 
hardness  causing  an  “unfurling’’  of  the  cilium,  as  discussed  later)  is  a  mechanism  for  reducing 
drag  during  the  cilium  return  stroke;  conversely,  the  full  extension  of  R^g  (increase  in  hardness) 
is  a  mechanism  for  maximizing  the  output  power  during  the  power  stroke.  The  x-acceleration  in 
figure  8  has  maxima  and  minima  at  the  boundaries  of  the  strokes.  The  use  of  Rayg  has  produced 
consistent  results.  Considering  the  area  about  the  zero  crossings,  the  thrust  efficiency  (the 
difference  between  the  positive  and  negative  areas  divided  by  the  total  area  under  the  curve)  is 
estimated  to  be  0.125. 
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Velocity*Length 


X-velocity  at  Ravg  *  Ravg  length  (in  y-z  plane) 
x  aq*  Biological  Cilia  - 17  um  length  -  19-Dec-201 1 


Figure  10.  Calculated  Variation  of  the  Variable  (Velocity  \  Ravg  in  pm2 Vs) 
of  the  Biological  Cilium  in  the  Revisited  TM  Data 


3.2  ROBOTIC  CILIUM 


Before  arriving  at  the  final  design  of  the  cilium  drive,  several  studies  were  carried 
out  to  better  understand  the  nature  of  cilium  motion.  These  included  a  computer-aided 
design  (CAD)  simulation,  the  development  of  a  hardware  D-cam  pendulum .  development 
of  an  early  drive  mechanism  using  orthogonal  bucket  handles,  and  a  robotic  cilium  design 
for  bending  motion  (see  appendix  sections  A. 3  and  A.4  for  details  of  these  studies). 


The  inset  in  figure  1 1  shows  a  set  of  independent  orthogonal  oscillatory  motions 
(IOOMs)  in  a  Cartesian  coordinate  system  that  has  been  used  to  generate  the  motion  in 
the  final  design  of  the  biorobotic  cilium.  It  has  two  roll  oscillations  in  orthogonal  x-  and 
z-axes.  Also,  a  twist  oscillation  in  they-axis  and  a  bending  oscillation  have  been  added 
to  account  for  the  variation  of  RavR  during  the  power  and  return  strokes. 


In  section  A.4  of  the  appendix  and  in  the  DVD  in  Bandyopadhyay  (2009),  it  is  noted 
that  these  IOOMs  have  the  ability  to  reproduce  the  motion  of  the  pectoral  fins  in  larger 
swimming  animals  as  well;  see  section  A.7.  Among  the  robotic  IOOMs,  the  x  and  z  roll 
motors  and  the  twist  motor  apply  moments,  and  the  push-pull  forces  are  applied  to 
generate  in-plane  bending. 


Figure  1 1  is  a  photograph  of  the  four-oscillating-motor  apparatus  producing  the 
actuations  shown  in  the  inset.  Each  motor  (Maxon)  has  position  encoders,  and  the  torque 
can  be  determined  from  the  input  current.  The  bending  oscillation  is  produced  by  a  pair 
of  sliding  flats  as  shown  in  figure  A-9  of  the  appendix;  a  flexible  tube  encapsulates  the 
push-pull  sliding  flats.  The  robotic  cilium  length  is  0.21  m. 

Figure  12  is  a  strobe  picture  sequence  of  cilium  position  during  the  return  stroke 
produced  in  the  apparatus  shown  in  figure  1 1 .  The  strobe  picture  may  be  compared  with 
figure  4b. 
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Figure  11.  Photograph  of  the  Scaled-Up  Model  of  a  Single  Robotic  Cilium 
(The  x-  and  z-motors  (roll)  lie  in  the  plane  of  the  horizontal  aluminum  frame  and  are  visible. 

The  y-motor  (twist)  and  the  bending  motor  He  underneath  the  flat  aluminum  frame.  Inset: 
Elementary  orthogonal  oscillations  of  a  cilium  protruding  from  a  surface  (x,  z).  Power  stroke 
lies  predominantly  in  the  (x,  y)  plane;  (y,  z)  is  the  transverse  plane.) 


Figure  12.  Strobe  Picture  of  the  Cilium  in  the  Apparatus  Shown  in  Figure  11 
(Return  stroke  is  shown;  -X  is  upstream  direction.) 
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3.3  COMPARISON  OF  ROBOTIC  CILIUM  AND  BIOLOGICAL  CILIUM 


In  this  section,  the  motions  of  the  scaled-up  robotic  cilium  are  compared  with  those  of  the 
biological  cilium.  Note  that  the  robotic  cilium  motion  in  air  is  being  evaluated  when  it  is  made 
to  reproduce  the  fluid-loaded  cilium.  The  beat  period  of  the  robotic  cilium  is  immaterial  and  is 
lower  by  a  factor  of  10  to  remove  inertia  effects.  The  units  in  the  biological  cilium  are:  length 
(pm),  velocity  (pm/s),  acceleration  (pm/s2),  curvature  (rad/pm),  and  torsion  (rad/pm).  The  units 
in  the  robotic  cilium  are  length  (m),  velocity  (m/s),  acceleration  (m/s2),  curvature  (rad/m),  and 
torsion  (rad/m). 


3.3.1  Comparison  of  Position 

In  figure  13a,  side  and  top  views  ({x(t),y(t)}  and  {x(t),  z(t) }  distributions,  respectively)  of 
cilium  position  during  the  beat  cycle  are  compared  for  the  biological  cilium  and  the  robotic 
cilium.  The  green  arrow  indicates  the  power  stroke,  and  the  red  arrow  indicates  the  return 
stroke.  The  patterns  in  the  biological  cilium  and  robotic  cilium  are  similar. 


Side  view 


Figure  13a.  Comparison  of  the  Biological  Cilium  (left  column)  and  Robotic 
Cilium  (right  column)  Positions — Side  View:  (x(t),  y(t)),  and  Top  View:  (x(t),  z(t)) 
(x-axis  is  in  the  horizontal  direction;  green  arrow  is  the  power  stroke;  red  arrow 
is  the  return  stroke.  Color  code:  Different  time  steps.  Biological  cilium  is  the 
UBC  version  ofTeunis  and  Mach  enter  (1994)  data.) 
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Figures  13b  and  13c  compare  the  biological  and  robotic  cilium  rotational  tracks;  in  figure 
13b.  20  positions  along  the  length  of  the  cilium  are  shown;  in  figure  13c.  3  locations  are  shown. 
The  topology  in  figure  13b  has  a  variable  precession  rate,  and  the  coupling  between  the  cilium 
beating  and  the  precession  rate  (Bouzarth  et  al..  2007)  that  the  fluid  particles  experience  has  a 
cyclic  variation.  The  topology  resembles  a  bent  horn,  but  is  not  easily  describable  by  a  complex 
Fibonacci  function  although  a  start  from  negative  integers  gives  the  bending. 

3D  Path  fits  with  respect  to  time 
Normalized  bio  and  mechanical  cilia  -  14-Dec-2011 


-0.2 


y  (vertical)  x  (propulsion) 

Figure  13b.  Comparison  of  Tracks  in  the  Biological  Cilium  (blue)  and  Robotic  Cilium  (red) 

3D  Path  fits  with  respect  to  time 
Normalized  bio  and  mechanical  cilia  -  14-Dec-2011 


Figure  13c.  Comparison  at  Cilia  Axial  Distances  of  S/L  =  0.33,  0.66,  and  1.0 
(The  base  (S/L  =  0)  is  indicated  by  the  dot;  blue:  biological  cilium;  red:  robotic  cilium.) 
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Figure  14  shows  the  projection  of  the  cilium  on  the  propulsion  plane.  (In  a  redesigned 
version  of  the  cilium.  the  bending  amplitude  has  been  increased,  which  would  allow  better 
agreement  with  the  biological  cilium  in  figure  14.)  The  cilium  is  precipitously  losing  its 
hardness  at  TS  >13.  It  is  suspected  that  in  the  biological  cilium.  the  material  property  is 
changing  abruptly  at  TS  >13.  over  the  beat  cycle.  The  torsional  pendulum  and  the  autonomic 
modeling  given  earlier  imply  that  the  cilium  is  acting  like  a  spring,  storing  and  releasing  energy 
during  the  return  and  power  strokes,  respectively.  This  view  will  be  synthesized  later. 


3D  plot  of  fitted  data 

Mechanical  Cilia  -  0  21  m  length  -  IB- Sep  2010 
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3D  plol  of  fitted  data 

Bio  Cilia  Data  17  am  length  -  10-Sep-2010 
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(a)  Robotic  Cilium 


(b)  Biological  Cilium  (Revisited  TM  Data) 


Figure  14.  Comparison  of  the  Robotic  and  Biological  Cilium  Positions  at 
Each  Time  Step  on  the  Longitudinal  Plane 

(The  red  dot  shows  the  base  location;  the  base  cilium  is  shifted  horizontally  to  avoid  clutter.) 


3.3.2  Comparison  of  Curvature  and  Torsion 

Define  r  (s,  t )  as  the  position  vector  of  a  point  on  the  surface  of  the  cilium,  w  ith  primes 
denoting  the  derivatives  with  respect  to  /,  and  x  denoting  the  cross  product;  1 1  as  magnitude, 
sign  is  retained.  The  following  equations  were  used  to  calculate  curvature  and  torsion  for  the 
biological  cilium  and  the  robotic  cilium  data  (Zwillinger,  1996).  Curvature  is  defined  as 


K'  = 


|r'X  /-"I 


Torsion  is  defined  as 

r'ir''Xr'”) 

T=~. - 7^- 

\r'Xr'  f 


(3) 


(4) 
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Curvature  measures  the  deviance  of  a  curv  e  from  being  a  straight  line  relative  to  the 
osculating  plane.  In  the  elementary  differential  geometry  of  curves  in  three  dimensions,  the 
torsion  of  a  curve  measures  how  sharply  it  is  twisting.  If  the  torsion  is  zero,  the  curve  lies 
completely  in  the  same  osculating  plane  (there  is  only  one  osculating  plane).  Note  that  the 
curvature  and  the  torsion  of  a  helix  are  constant;  when  they  are  not  constant,  the  geometry  is  not 
helical. 

Curvature  K  and  torsion  T  in  the  robotic  and  biological  cilium  are  compared  in  figures  15 
and  16.  The  biological  ciliunvs  curvature  and  torsion  shown  are  from  the  revisited  TM  data. 

The  white  lines  in  figures  1 5  and  1 6  indicate  the  boundaries  of  the  power  and  return  strokes;  the 
return  stroke  lies  between  the  white  lines.  These  results  are  different  from  those  of  Teunis  and 
Machemer  (1994).  In  our  definition,  the  power  and  return  strokes  are  of  equal  extent — not  1/7  to 
2/7  for  power,  as  given  by  Teunis  and  Machemer.  In  figures  1 5  and  16,  TS  1  and  42  are 
assigned  phase  of  0  and  2 n.  that  is,  the  phase  at  the  /7th  time  step  is  <}>  =  Inn  /  N ,  where  N(=4\) 
is  the  total  number  of  time  steps. 

In  figures  1 5  and  1 6,  for  the  robotic  cilium  data  after  lowpass  filtering,  0  <  K  (rad/m)  < 

1 0.8 1 1 7,  -1 .78  <  T (rad/m)  <  35.87.  For  the  biological  cilium  data  after  lowpass  filtering,  0  <  K 
(rad/pm)  <  0.21 78,  -0.0729  <  T (rad/pm)  <  0.4834.  The  minimum  values  are  in  blue  and  the 
maximum  values  are  in  red.  The  biological  cilium  and  the  robotic  cilium  have  approximately 
similar  curvature  and  torsion  distributions. 

The  calculation  of  curvature  is  consistent  w  ith  the  conclusion  of  Teunis  and  Machemer  in 
that  there  is  little  curvature  during  the  power  stroke.  This  finding  can  be  further  verified  by 
looking  at  the  position  data  in  figure  14.  The  torsion  calculation  also  matches  the  conclusions 
above.  In  the  biological  cilium,  curvature  is  high  near  the  base  after  the  return  stroke  has  started. 
Torsion  reaches  a  peak  approximately  halfway  along  the  length  when  the  power  stroke  ends 
(marked  as  point  P  in  figure  16b). 

Curvature  and  torsion  limits  widen  if  the  cilium  length  is  reduced  from  50  to  17  pm.  Teunis 
and  Machemer  (1994)  state  that  curvature  and  torsion  are  small  during  the  power  stroke  (that  is 
the  cilium  is  stiff).  They  state  that,  “Both  curvature  and  torsion  return  to  minimal  values  by  the 
beginning  of  the  power  stroke”  (see  figures  1 5  and  1 6). 
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Figure  15.  Comparison  of  the  Curvature  of  the  Robotic  Cilium  (a)  and  the 
Biological  Cilium  (b)  Using  the  Revisited  TM  Data 
(See  text  for  color  code.) 
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Figure  16.  Comparison  of  Torsion  in  the  Robotic  Cilium  (a) 
and  the  Biological  Cilium  (b) 

(See  text  for  color  code;  P  is  a  point  of  inflection 
on  the  biological  cilium  (see  figure  1 7).) 


4.  SYNTHESIS  OF  RESULTS 


4.1  MODELING  THE  INCREASE  AND  DECREASE  OF  THE  BIOLOGICAL 

CILIUM’S  HARDNESS 

I'he  curvature  and  torsion  results  are  synthesized  in  figure  17  to  understand  the  behavior  at  the 
transition  from  the  power  to  the  return  stroke  when  the  hardness  of  the  cilium  changes  precipitously. 

The  model  schematic  in  figure  17  shows  projections  of  the  cilium  onto  the  propulsion  plane 
at  TS  9,  14.  and  16.  These  are  copies  from  figure  14b.  Between  TS  13  and  14,  the  time 
derivative  of  the  angle  subtended  by  the  cilium  at  the  base  changes  sign  (the  cilium  positions 
converge  in  figure  14b.  This  is  depicted  in  figure  17  by  the  change  in  the  sign  of  the  moments 
applied  at  the  base  m.  The  directions  of  velocities  V of  the  cilium  and  of  the  drags  D  acting  on  it 
are  shown.  At  TS  14,  while  the  lower  part  of  the  cilium  turns  anti-clockwise,  the  distal  part  turns 
in  the  opposite  direction.  As  a  result,  a  point  of  inflection  is  produced  in  between  (marked  as  P) 
where  the  finite  radius  of  curvature  r  changes  sign.  The  cilium  precipitously  droops  after  TS  14, 
which  indicates  a  decrease  in  hardness.  The  cilium  is  straight  during  the  power  stroke. 

Therefore,  the  question  is:  How  does  the  hardness  drop  and  recover  in  every  beat  cycle?  This 
process  is  modeled  in  the  insets  P.  PI.  and  P2  in  figure  17a  and  in  figure  16b  by  applying 
fracture  mechanics. 

Figure  16b  shows  that  at  point  P  curvature  is  not  large,  but  torsion  reaches  a  peak  because  the 
cilium  is  being  twisted  in  the  transverse  plane.  This  is  shown  in  the  insets  in  figure  17a.  Inset  P 
shows  that  at  P.  the  neutral  plane  of  stress  undergoes  a  reversal  in  compressive  C  and  tensile  T 
stresses.  The  two  ends  of  P  also  experience  a  reversal  in  the  sign  of  the  applied  torque  r  between 
the  cilium  base  and  the  distal  point.  Inset  PI  shows  that,  as  a  result,  two  otherwise  parallel  circular 
cross-sections  shear  and  open  ajar  due  to  the  switching  of  the  regions  labeled  (’and  T. 

The  literature  suggests  that  in  the  9+2  axoneme  the  central  microtubule  pair  is  responsible 
for  hardness  control.  It  can  also  be  surmised  that  motor  protein  molecules  climb  up  the 
microtubule,  find  attachment  spots,  and  lock  (Cordova  ct  al..  1992).  This  process  is  modeled  in 
figure  1 7b  and  inset  P2  in  figure  1 7a.  The  circular  part  of  inset  P2  shows  the  central  microtubule 
pair  ( CMau  and  CA//,/^)  and  the  cross-bridge  link  (CBpu)-  The  presence  of  many  such  cross¬ 
bridge  links  along  the  central  microtubule  pair  increases  the  moment  of  inertia  /,  increasing  the 
property  GI.  which  is  responsible  for  resistance  to  torsion,  where  G  is  the  modulus  of  torsional 
rigidity  of  the  material.  (The  cross-bridges  are  similar  to  cross-links  in  polymers  where  they  are 
closer  when  unstressed,  and  the  polymer  chain  is  straightened  when  the  links  are  pulled  apart  by 
the  application  of  stress.) 

It  is  modeled  that,  when  torsion  at  P  reaches  the  maximum  value,  the  cross-bridge 
attachment  cracks;  this  is  shown  in  the  boxed  part  of  inset  P2.  The  cross-bridge  at  location  P  at 
TS  14  is  the  proverbial  “last  straw” — when  the  cross-bridge  breaks,  there  is  no  hardness  left  for 
the  power  stroke  to  continue.  Inset  P2  shows  how  a  crack  develops  (the  sphere  is  the  foot  of  the 
cross-bridge)  at  the  sites  where  there  the  cross-bridge  is  attached.  The  cross-bridge  motor 
protein  and  the  host  microtubule  site  are  of  dissimilar  materials  conformationally  held  in  place. 


27 


on  the  Cilium  at  the  Beat  Phase  of  Neutral  Equilibrium  (See  Inset  P) 

(Symbols:  9,  14  and  16  are  time  steps;  m  and  -m  are  moments  in  the  jc-j  plane  at  the  base;  D, 
V,  and  r  are  drag,  velocity,  and  the  radius  of  curvature;  the  prime  sign  (')  denotes  values  near 
the  base;  C  and  T  are  compressive  and  tensile  stresses;  and  N  =  neutral  plane.  The  cilium  is 
shifted  along  X  with  time  step  to  avoid  clutter.  Left  inset:  CMai4  is  the  central  microtubule  a  at 
TS 14;  CMbi4  is  the  other  centra!  microtubule  of  the  pair  at  TS  14;  CBp)4  is  the  cross-bridge  at 
the  point  of  inflection  P  at  TS  14.  The  cross-bridge  link  can  be  expected  to  fail first  at  the 

point  marked  by  the  red  dot  in  inset  P.) 


~T 


r~D 


Figure  1 7b.  Modeling  of  Cross-Bridge  Link  Between  the  Pair 
of  Central  Microtubules  (the  2  of  the  9  +  2) 

(Only  the  cross-bridge  at  point  P  is  shown;  this  is  the  only  one  that  cracks  because  torsion 
reaches  high  values  at  this  point;  other  cross-bridges  add  to  the  hardness,  but  they  do  not 
crack  because  the  torsion  is  not  high  (the  region  near  the  cilium  tip  is  ignored).  Broken  lines 
represent  later  times.  Curvature  is  small  near  point  P  at  TS  14;  therefore,  planar  vibrations 

have  a  negligible  effect  on  cracking.) 
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One  expects  cracks  to  develop  where  the  cross-bridge  is  attached  to  the  microtubule  and  not  at  its 
apex  because  the  microtubule  attachment  site  has  to  be  vacated  fully  for  the  new  beat  cycle  to  resume 
accepting  a  new  cross-bridge  attachment.  Our  modeling  considers  the  cross-bridge  cracking  at  the 
attachment  site  at  TS  14  at  point  P.  The  model  includes  cross-bridges  all  along  the  length  of  the 
microtubules,  but  only  the  one  at  the  point  P  at  TS  14  cracks  at  the  contact  site  and  needs  replacing. 


4.2  MODELING  OF  THE  CROSS-BRIDGE  DETACHMENT/SOFTENING 

As  per  Karp  (2005),  one  end  of  the  protein  attaches  to  one  microtubule  and  the  other  end  to 
another  microtubule,  creating  a  cross-bridge  structure.  Cordova  et  al.  (1992)  give  an  electrostatic 
model  of  the  cross-bridge  attachment  and  detachment  although  not  in  the  context  of  the 
paramecium's  cilium  beating.  Here,  fracture  mechanics  are  applied  to  model  the  detachment. 

The  cross-bridge-microtubule  contact  is  modeled  to  be  brittle  (like  glass),  rather  than  plastic  (like 
metal);  in  other  words.  Griffith's  law  (1920)  applies  and  Irwin's  (1948)  does  not.  Griffith's  law 
states  that  the  product  of  fracture  stress  and  the  square  root  of  the  crack  length  remains  constant 
as  the  crack  progresses  until  the  difference  between  the  surface  and  elastic  energy  reaches  the 
peak  value,  at  which  point  failure  occurs. 

Griffith's  relationship  may  be  applied  as  follows: 


°>  =~r- 

y]a 

(5) 

■=-t- 

(6) 

where  Of  is  the  stress  at  fracture,  a  is  the  crack  length,  C  is  a  constant,  E  is  Young's  modulus  of 
elasticity  (=  2  GPa  for  polyethylene;  =  (10  to  30)  x  10y  N/m2  for  myosin  (Adamovic  et  al., 
2008)),  and  y  is  the  surface  energy  density  (=  0.72  J/m  for  water,  1 .0  J/m“  for  glass,  and  >2 
J/irf  for  polymers;  bulk  cilium  is  modeled  as  a  water-like  polymer).  Other  parameters  are: 
equilibrium  spacing  of  the  cross-bridge  =  10  nm.  and  the  radius  of  the  actin  filament  =  5.5  nm 
(Cordova  et  al.,  1992).  Assume  the  value  of  a  to  be  ( 5.5;r )  nm  (hemispherical  arc  length),  E  (= 
0.50  GPa)  of  human  tendon  (Lauder  et  al.,  201 1)  to  be  applicable,  ( A ,  the  hemispherical  area 
2;r(5.5  nm)2;  and  note  that  1  J  =  1  Nm;  1  Pa  =  1  N/m2;  G  =  109;  1  n  =  10'9;  1  p  =  10'12;  y  =  10'24, 
then  fracture  stress  (aj)  would  be 

I  2(0.5  xlO9) 

V  ^25.5  x  10' * 

=  V2  x  1016 

=1.36  x  108  N/m2. 

The  force  F  required  to  fracture  one  of  the  foot  anchors  of  the  cross-bridge  =  2.6  x  1 0‘8  N  =  26 
nN  =  26  x  103  pN. 
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Earth's  gravity  exerts  a  force  of  10  -  100  pN  on  cells  in  bio-physical  systems.  In  their 
artificial  gravity  experiments  with  paramecia.  Guevorkian  and  Vales  (2006)  show  that  they 
respond  to  gravity  and  regulate  swim  speed  primarily  due  to  the  buoyancy  of  the  cell.  In  their 
thermal  cum  electrostatic  modeling  of  the  motion  of  a  single  motor  protein  (myosin)  molecule. 
Cordova  et  al.  (1992)  assume  a  maximum  cross-bridge  force  of  5  pN.  So.  a  cross-bridge 
detachment  force  of  26  x  104  pN  is  too  high,  and  5  pN  is  more  reasonable. 

On  the  other  hand.  Cordova  et  al.  (1992)  state  that  (1 )  the  attachment  and  detachment  sites 
of  the  cross-bridge  differ  by  the  presence  of  ATP  when  detached;  and  that  (2)  the  “binding  of 
ATP  weakens  the  attachment  of  the  cross-bridge  to  the  fiber.”  (ATP  (adenosine-tri-phosphate) 
transports  chemical  energy  within  cells  for  metabolism:  it  is  recycled;  it  is  also  a  messenger 
molecule;  it  is  sometimes  called  "the  molecular  unit  of  currency.”) 

In  linear  elastic  fracture  mechanics.  Griffith's  theory  (equations  (5)  and  (6))  states  that  free 
energy  is  given  by  the  difference  between  the  surface  energy  and  the  elastic  energy,  and 
increases  with  the  crack  length.  Failure  occurs  when  the  free  energy  attains  a  maximum  value  at 
a  critical  crack  length  beyond  which  the  free  energy  decreases  by  increasing  the  crack  length. 

The  free  energy  in  Griffith's  theory  may  be  corrected  due  to  the  hydrolysis  of  ATP  at  the 
detachment  site  when  the  chemical  energy  is  released.  One  can  back-calculate  what  this  should 
be  as  follows; 

The  force  Fthat  leads  the  cross-bridge  to  detach  =  [A.  the  hemispherical  attachment  area  x  aj\ 

=  2n(5.5  x  10’9)2  a,  (N)  =  190  x  10l8<r/(N).  For  this  force  F  to  be  5  pN,  <r,  =  5  x  lO12/ (1 .876 
x  10'15)  =  2.665  x  103  N/m2  =  2.665  kPa.  Assuming  as  before  E  =  0.5  x  109  N/m\  the  surface 
energy  density  y  would  have  to  be  (1  J  =  1  Nm)  3.856  x  1 0'1"  J/m2=  38.56  nJ/m  . 

The  ATP  hydrolysis  has  to  convert  the  chemical  energy  at  the  detachment  site  to  make  the 
surface  energy  density  to  be  38.56  nJ/nV  for  the  force  F required  (produced  practically  at  a  point) 
to  detach  the  cross-bridge  to  be  5  pN.  The  total  energy  (=  surface  energy  -  strain  energy  = 

2yuL  -  <j , '  m'  L/(2E)  =  yuL .  where  L  is  the  third  dimension  of  crack,  and  strain  energy  is  the 

elastic  energy  released  in  the  fracture  area  due  to  the  fracture;  the  strain  energy  comes  from  the 
torsion)  spent  per  detachment  is  75  x  10  24  J  =  75  yJ  (yocto  joules). 

Attributing  the  total  energy  to  the  chemical  energy,  the  chemical  energy  spent  is  estimated 
as  follows: 

One  mol  has  6.02  x  10  ‘  molecules;  the  energy  released  per  ATP  mol  is  45.6  U.  Therefore, 
the  energy  released  per  molecule  of  ATP  is  76  x  103  yJ,  which  can  power  1000  cross-bridge 
detachments.  For  3000  cilia  in  one  paramecium,  3  ATP  molecules  are  needed  per  beat 
cycle  to  detach  one  cross-bridge  in  each  cilium;  so,  at  14  Hz.  42  ATP  molecules  are  needed 
per  second  (fewer  are  needed  if  the  absolute  viscosity  drops). 

A  full  modeling  of  the  ATP  hydrolysis  process  is  needed  to  quantitatively  complete  the 
evaluation  of  the  mechanism.  The  hardness  control  model  should  include  the  mechanisms  of 
softening  at  the  detachment  point  due  to  the  presence  of  ATP,  subsequent  enhanced  bonding  due 
to  ATP  hydrolysis,  the  effect  of  torsion,  and  sliding  between  the  microtubules. 
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5.  DISCUSSION 


5.1  DEVELOPMENT  OF  THREE-DIMENSIONAL  DEFORMATION  IN  CILIUM 

AND  THE  ROLE  OF  CRITICAL  HARDNESS 

The  role  of  the  modulus  of  torsional  rigidity  in  the  generation  of  torsion  (=  three- 
dimensionality)  in  the  cilium  was  investigated  experimentally  and  is  presented  in  section  A. 5  of 
the  appendix.  Three-dimensionality  can  lead  to  stability  of  cruising  and  facilitation  of  food 
gathering.  The  paramecium  cruises  in  a  stable  corkscrew  path  about  an  imaginary  axis  while 
rotating  about  its  own  body's  axis.  Compared  to  two-dimensional  motion,  the  three-dimensional 
motion  of  the  cilium  imparts  more  momentum  to  the  water  and  allows  the  thrust  vector  to  be 
slanted  to  the  plane  of  the  power  stroke  (  figure  7).  In  three-dimensional  motion,  all  phases  of  the 
motion  contribute  to  the  production  of  an  inw  ard  acceleration  of  the  tiny  tornadoes.  These 
tornadoes  (or  their  agglomerations)  could  allow  the  food  elements  to  be  drawn  in  as  well. 

Section  A. 5  of  the  appendix  show  s  the  curvature  and  torsion  in  several  models  of  long  cilia 
made  of  materials  of  different  hardness  but  operated  by  the  same  type  of  actuator,  consisting  of 
three  doublets  (pairs  of  push-pull  strings)  spaced  120°  apart.  The  hardness  varied  from  near  zero 
to  a  Shore  hardness  of  A53-63.  When  hardness  was  very'  low,  a  helix  was  produced  (a  helix  has 
a  constant  curv  ature  and  torsion).  In  contrast,  in  the  high-hardness  samples,  only  planar 
curvature  was  produced  and  torsion  was  zero  (there  was  no  off-plane  curvature).  The  conclusion 
was  that  Shore  hardness  should  be  below  a  certain  critical  value  for  torsion  to  be  manifested.  At 
the  critical  value  of  hardness,  the  cilium  would  he  sensitive  to  any  attempt  at  reduction  or 
increase  in  hardness.  The  detachment  of  only  one  cross-bridge  in  the  biological  cilium  right 
after  TS  14  at  a  strategic  location  on  the  cilium  (location  P),  and  the  attachment  of  a  new  one  at 
the  same  location  P  in  the  follow  ing  phases  would  seem  to  allow  control  of  the  hardness  at  a 
critical  level  for  low  cost  (a  behavior  akin  to  a  high  quality  factor  in  resonant  systems). 


5.2  BIO-PHYSICAL  SELF-REGULATION  OF  CILIUM  HARDNESS 

The  dynamic  systems  analysis,  the  modeling  of  hardness  control,  and  the  preceding 
discussion  on  the  value  of  designing  a  cilium  of  critical  hardness  point  to  the  existence  of  the 
phase  of  neutral  equilibrium  and  minimal  energy  expenditure.  The  authors  deem  this  to  be  self¬ 
regulation. 

It  is  known  that  when  a  paramecium  hits  an  obstacle,  it  backs  away  at  an  angle  and  proceeds 
in  a  new  direction  (dePeyer  and  Machemer.  1978).  The  collision  with  the  obstacle  can  be 
modeled  via  the  amplitude  and  duration  of  the  extra-cellular  impulse  term  lex,  in  equation  (1 ). 
The  cilium  motion  follows  a  limit  cycle.  The  return  path  to  the  limit  cycle  is  unique  for  a  given 
value  of  Iex,  and  the  phase  of  the  system  when  Iex,  is  applied. 
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The  3000  cilia  in  a  paramecium  are  each  a  nonlinear  oscillator.  The  self-referential  phase 
reset  (SPR)  property  of  equation  (1)  allows  for  synchronization  between  actuators  that  receive 
the  extracellular  impulse  (Kazantsev  et  al.,  2003  and  2004:  and  Bandyopadhyay  et  al.,  2008).  So 
far.  the  evidence  for  this  process  is  conflicting. 

Yeh  et  al.  (1990)  have  observed  spectral  peaks  at  370  kHz  (attributed  to  the  myosin  thick 
filament)  and  a  relaxation  time  of  5  ps  (attributed  to  cross-bridge  motion).  To  straighten  the 
cilium  during  the  return  stroke,  the  paramecium  fuel  cell  needs  to  provide  energy.  The  power 
stroke  makes  use  of  this  energy.  At  the  resonant  frequency,  maximum  current  flows  through  R 
in  an  RLC  circuit.  If  C is  10  nF  and  L  is  10  nH,  the  resonance  frequency  (f,  =  (\I(2tc))  * 

[1  N(LC)])  is  1 6  MHz.  This  is  x  10h  of  14  Hz,  the  beat  frequency  of  the  cilium  in  water  at  20°  C. 
If  more  cross-bridges  link  the  microtubule  pair  increasing  C  and  L  to  10  pF,  and  10  pH, 
respectively,  the  resonant  frequency  would  be  1 6  kHz,  which  is  x  1  O'*  of  14  Hz — a  value  much 
closer  to  that  observed  by  Yeh  et  al.  Parametric  oscillators  with  second-order  nonlinear 
interactions  can  produce  lower  frequencies  from  two  higher  frequencies  by  differencing  their 
frequencies:  however,  there  is  no  known  evidence  that  this  is  happening. 

If  <f>  is  torsional  angle,  the  equation  of  torsional  motion,  as  per  a  spring-mass-damper  model 
is 


10  +  2^a>n0  +  (Tk)<j>  =  0 , 


where 


is  the  natural  frequency  of  torsional  oscillation  of  the  system  (  a>0  =  2 nfa ), 


and  Q  =  — is  the  damping  ratio,  Tc  is  the  torsional  damping  (apart  from  viscous  damping, 
2yjCk)I 

there  may  be  structural  damping  also,  which  is  hysteretic),  T  k  is  the  torsional  spring  constant, 
(applied  torque  ( r k</>)  +  damping  torque  ( r c<f>)  is  the  total  torque),  and  /  is  the  moment  of  inertia. 


The  description  of  a  novel  hemispherical  motor  (‘‘frictionless")  cilium  actuator  is  given  in 
section  A.4  of  the  appendix.  The  present  authors  have  successfully  operated  the  IOOMs  of  the 
robotic  cilium  w  ith  the  real-time  analog  solution  of  equation  (1),  as  described  in  Bandyopadhyay 
et  al.  (2008  ).  This  integrated  operation  of  the  controller  and  the  actuator  is  discussed  in  section 
A. 6  of  the  appendix.  These  two  developments  take  us  one  step  closer  to  fuller  robotic  simulation 
of  cilium  motion. 
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5.3  RELATIONSHIP  WITH  LARGER  SCALE  NATURAL  SWIMMING:  A 
SYNTHETIC  VIEW 

As  per  Llinas  (2001),  Llinas  and  Yarom  (1081),  and  Llinas  et  al.  (2004): 

1.  Motion  control  neurons  follow  self-regulating,  dynamic  systems  principles. 

2.  The  neurons  can  be  clustered  together  on  demand  (self-referential  phase  reset  (SPR) 
property  (Kazantsev  et  al.,  2004). 

3.  Neuronal  development  in  animals  parallels  that  of  biomechanics,  the  former  being  mature 
in  motion  control  over  many  species  and  remaining  unchanged. 

Extending  these  principles  further,  if  sensors,  actuators,  and  their  control  follow  the  same 
dynamic  systems  principles  to  remain  in  persistent  synchrony  with  the  environment,  then  they 
home  faster  on  moving  targets  if  they  are  assumed  to  have  a  10-15%  stereoscopic  bias  in  sensing 
and  actuation  (Bandyopadhyay,  2012a). 

A  flapping  fin  is  lift  based  and  it  flaps  cross  stream  to  the  thrust.  On  the  other  hand,  a  cilium 
is  drag  based  and  it  flaps  both  in  the  plane  of  thrust  and  in  the  transverse  plane.  How  does  the 
transition  from  small  to  large  and  from  drag-based  to  lift-based  propulsion  take  place? 

Ellington  (1999)  has  considered  the  wing  flapping  in  insects  to  be  a  resonant  oscillation  that 
has  low  friction  and  high  amplitude.  Triantafyllou  and  Triantafyllou  (1995)  have  shown  that 
varieties  of  fish  flap  their  caudal  fins  in  a  range  of  preferred  frequency  called  the  Strouhal 
number  (0.2  to  0.4).  Note  that  the  length  scale  in  Strouhal  number  is  given  by  the  amplitude  of 
oscillation,  which  can  be  expected  to  reach  a  large  value  in  resonant  oscillation.  In  recent  work, 
the  present  authors  have  been  able  to  operate  penguin-wing-like  flapping  fins  in  large  scale  for 
ocean  engineering  purposes  while  flapping  at  the  above  range  of  Strouhal  number 
(Bandyopadhyay  et  al.,  2011). 

Paramecia  have  pressure-sensitive  mechanoreceptors  (dePeyer  and  Machemer.  1978  and 
1983),  sliding  microtubules  that  produce  curvature  in  cilia,  phased  beating  of  cilia  along  their 
body  surface  that  produce  the  animal's  tumbling  motion  in  a  corkscrew  path,  and  a  9+2 
axenome  architecture.  Fish  have  pressure-sensitive  lateral  lines,  hemitrich  bone  pairs  in  their 
pectoral  fins  for  producing  curved  shapes  (Lauder  et  al..  2011),  central  pattern  generators  for 
producing  tail  beat,  and  1 1  hemitrich  bones  on  their  pectoral  fin  spaced  with  a  common  origin 
(Lauder  and  Drucker,  2004).  These  similarities  were  assumed  to  owe  their  origin  to  co-evolution 
of  sensors,  actuators,  and  control,  and  a  common  oscillatory  framework  of  a  dynamic  system.  In 
section  A. 7  of  the  appendix  (see  Bandyopadhyay,  2009),  the  outline  of  a  notional  synthetic 
description  vividly  shows  the  common  lineage  of  large  and  small  swimming  animals. 

Each  of  the  tiny  individual  cilium  propulsors  appears  to  have  the  intrinsic,  dormant 
kinematic  and  structural  building  blocks  required  to  optimize  at  both  low  and  high  Reynolds 
numbers.  The  tentative  conclusion  is  that  the  integrated  autonomous  sensory,  actuating,  and 
controlling  architecture  of  the  paramecium  is  invariant  of  Reynolds  number. 
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6.  CONCLUSIONS 


Modeling  and  analysis  of  the  motion  of  the  paramecium's  cilium  have  been  carried  out,  with 
a  focus  on  the  motion's  oscillatory  nature.  Theoretical  modeling  included  considerations  of  a 
torsional  pendulum,  dynamic  systems,  and  fracture  mechanics.  A  revised  analysis  of  Teunis  and 
Machemer's  (1994)  photo-microscopy  of  biological  cilium  has  been  compared  with  a  robotic 
simulation  of  the  biological  cilium.  The  role  of  torsional  rigidity  in  three-dimensional 
deformation  has  been  examined.  The  aspects  of  integrating  sensors,  controller,  and  actuators  that 
remain  invariant  with  Reynolds  number  have  been  robotically  simulated  and  discussed. 

The  following  conclusions  can  be  drawn: 

1 .  The  variation  of  the  cilium  beat  frequency  with  fluid  viscosity  can  be  modeled  by 
describing  the  beat  frequency  as  the  natural  frequency  of  a  torsional  pendulum. 

2.  The  cilium  track  is  a  limit  cycle  describable  by  dynamic  systems  theory.  There  exists  a 
phase  of  neutral  equilibrium  lying  at  the  boundary  where  the  power  stroke  changes  to  the  return 
stroke  and  the  cilium  velocity  is  zero. 

3.  Several  long,  mechanical  cilia  have  been  fabricated,  and  it  was  found  that  their  hardness 
has  to  be  below  a  certain  level  for  torsion  to  be  manifested. 

4.  The  three-dimensional  motion  of  each  cilium  produces  a  tiny  tornado  of  fluid  having  an 
inward  acceleration.  The  inward  acceleration  is  produced  during  both  the  power  and  return 
strokes. 

5.  I'he  local  angle  that  the  cilium  orbital  velocity  and  the  cilium  position  makes  with  the 
propulsion  axis  remain  nearly  constant,  along  the  length  of  the  cilium  during  the  power  stroke, 
but  it  varies  widely  during  the  return  stroke. 

6.  The  average  radial  length  of  the  cilium  in  the  transverse  plane  remains  large  during  the 
first  two-thirds  of  the  power  stroke,  dropping  precipitously  during  the  remaining  one-third.  The 
average  cilium  radial  length  increases  during  the  return  stroke.  This  variation  of  the  average 
radial  length  should  help  increase  thrust  during  the  power  stroke  and  lower  drag  during  the  return 
stroke. 

7.  The  zero  crossings  of  the  product  of  the  average  radial  length  of  the  cilium  in  the 
transverse  plane  and  of  the  cilium  velocity  are  used  to  define  the  extents  of  the  power  and  return 
strokes.  The  power  stroke  is  0.50  of  the  beat  period,  which  is  larger  than  the  0.14  to  0.28  range 
given  by  Teunis  and  Machemer  (1994). 

8.  Using  four  independent  orthogonal  oscillatory  motor  drives,  a  robotic  cilium  has  been 
constructed  whose  curvature  and  torsion  distributions  are  in  fair  agreement  with  those  of  the 
biological  cilium.  The  robotic  cilium  has  also  been  controlled  by  analog  olivo-cerebellar 
dynamics. 
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9.  A  bio-physical  model  of  self-regulation  of  cilium  hardness  is  given  where  cross-bridge 
links  between  the  microtubule  pairs,  ATP  molecules,  and  torsion  play  a  key  role  in  hardness 
control.  Fracture  mechanics  has  been  applied  to  the  cross-bridge  detachment  at  the  phase  of 
neutral  equilibrium  at  the  cilium  location  where  a  point  of  inflection  is  created.  Synthesis  of  the 
results  suggests  that  the  cilium  is  torsionally  wound  (like  a  spring)  during  the  return  stroke  and 
the  energy  is  expended  during  the  power  stroke. 

10.  The  similarities  in  the  sensory,  actuating,  and  control  architecture  of  paramecia  and  fish 
are  noted. 
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7.  FUTURE  WORK 


Future  work  should  he  carried  out  in  the  following  areas: 

1.  Polymer  Nano-Micrometer  Scale  Mechanical  Cilium:  Such  a  cilium  can  be  fabricated 
as  follows.  It  is  knowTi  that  when  solvated  polymer  is  ejected  through  a  positively  charged 
conducting  syringe  by  a  pump,  the  laminar  jet  flow,  when  allowed  to  dry,  forms  a  cylindrical 
fiber.  The  instability  of  such  jets  is  considered  by  Rallison  et  al.  (1995).  Hohman  et  al.  (2001), 
and  Shin  et  al.  (2001).  From  a  w  ide  range  of  polymeric  materials,  such  an  electrospinning 
technique  can  be  used  to  make  cylindrical  nanofibers  of  length  1 7  pm  with  a  diameter  of  less 
than  1  pm  (see  Lu  et  al.  (2009),  Burger  et  al.  (2006),  Shin  et  al.  (2001),  Reneker  and  Chun 
(1996),  Reneker  et  al.  (2007),  Ramaswamy  et  al.  (201 1),  Ojha  et  al.  (2008),  and  McCullen  et  al. 
(2007)).  Such  cilia  would  have  the  Reynolds  number  of  the  biological  cilium.  The  mechanical 
rendition  consists  of  a  pair  of  such  electrospun  polymer  fibers,  closely  spaced  but  attached  at  the 
distal  point.  The  base  of  the  assembly  is  mounted  on  a  pair  of  stacked  MEMS  (micro-electro¬ 
mechanical  systems)  tables.  The  table  assembly  is  oscillated  orthogonally  as  in  figure  1 1.  The 
pair  of  nanofibers  is  given  a  push-pull  oscillation  to  generate  bending.  Using  another  MEMS 
table  underneath  the  stacked  pair,  the  entire  assembly  described  above  is  given  a  twisting 
oscillation  as  in  figure  1 1  (Bandyopadhyay.  20012b). 

2.  “Frictionless  Actuation”  with  a  Hemispherical  Motor  Drive:  Using  programmable 
logic  circuits  and  the  0.21-m-long  robotic  cilium  (or  the  polymer  nano-cilium)  base  torque 
actuation  with  the  hemispherical  motor  drive  can  be  explored.  Such  a  gearless,  noncontact  drive 
can  be  expected  to  be  nearly  frictionless.  This  would  simulate  the  condition  of  resonant 
oscillation. 

3.  Disturbance  Rejection  Verification:  Further  verification  of  the  limit  cycle  nature  of  the 
cilium  trajectory  can  be  made  using  the  hemispherical  motor.  Disturbances  can  be  applied  to  the 
cilium  to  see  the  trajectory  it  takes  to  return  to  the  ECO.  The  return  paths  can  be  compared  with 
the  theoretical  trajectories.  This  experiment  would  verify  the  autonomic  nature  of  cilium 
swimming. 

4.  Three-Dimensional  Hydrodynamic  Modeling:  Such  modeling  would  bring  us  one  step 
closer  to  exploring  three-dimensional  metachronism. 

5.  Origin  of  Metachronism:  Application  of  dynamic  systems  theory  is  a  new  approach  to 
the  understanding  of  cilium  autonomy.  A  parallel  exploration  of  metachronism  due  to 
hydrodynamic  coupling  of  neighboring  cilia  and  that  due  to  term  Iext  in  equation  (1)  could  be 
carried  out. 

6.  Modeling  Torsional  Effects:  The  effect  of  torsion  on  the  cross-bridge  linking  between 
microtubules  could  be  theoretically  modeled  to  understand  how  hardness  can  be  varied  to  store 
and  release  mechanical  energy  in  a  torsional  oscillatory  system  with  minimal  energy 
expenditure. 


37  (38  blank) 


8.  BIBLIOGRAPHY 


Adamovic,  I.,  S.  M.  Mijailovich.  and  M.  Karplus  (2008),  “The  Elastic  Properties  of  the 

Structurally  Characterized  Myosin  II  S2  Subdomain:  A  Molecular  Dynamics  and  Normal 
Mode  Analysis,’*  Biophysics  Journal,  vol.  04,  pp.  3770-3789. 

Au,  W.  (1993),  The  Sonar  of  Dolphins,  Springer-Verlag.  New  York. 

Bandyopadhyay,  P.  R.  (2012a),  “Handedness  Helps  Homing:  Smart  Materials,  Structures  and 
Systems,”  CIMTEC  2012,  June  10-14,  2012.  Montecatini  Terme,  Italy. 

Bandyopadhyay.  P.  R.  (2012b),  “A  Nano-Micrometer  Scale  Mechanical  Cilium  for  Dynamic 
Transduction,”  Invention  Disclosure.  Navy  Case  (to  be  submitted).  Naval  Undersea 
Warfare  Center  Division,  Newport.  RI. 

Bandyopadhyay,  P.  R.  (201 1),  “An  Acousto-Optical  Method  of  Encoding  and  Visualization  of 
Underwater  Space.”  Invention  Disclosure.  Navy  Case  Number  101454,  submitted 
November  201 1  (patent  pending). 

Bandyopadhyay.  P.  R.  (2009).  “Origin  of  Propulsion:  A  Synthetic  View,”  DVD  (edited  with 
audio).  Naval  Undersea  Warfare  Center  Division,  Newport,  Rl  02841. 

Bandyopadhyay,  P.  R.  (2005),  “Trends  in  Biorobotic  Autonomous  Undersea  Vehicles,”  IEEE 
Journal  of  Oceanic  Engineering,  vol.  30.  no.  1.  pp.  109-139. 

Bandyopadhyay,  P.  R.  (2004).  “Biology-Inspired  Science  and  Technology  for  Autonomous 

Underwater  Vehicles,”  IEEE  Journal  of  Oceanic  Engineering,  vol.  29,  no.  3.  pp.  542-806. 

Bandyopadhyay.  P.  R.,  H.  A.  Leinhos.  J.  D.  Hrubes.  N.  Toplosky.  and  J.  Hansen  (201 1), 
“Turning  of  a  Short  Length  Cable  Using  Flapping  Fin  Propulsion.”  IEEE  Journal  of 
Oceanic  Engineering,  vol.  36,  pp.  571-585. 

Bandyopadhyay,  P.  R.,  S.  N.  Singh,  D.  P.  Thivierge.  A.  M.  Annaswamv.  H.  A.  Leinhos,  A.  R. 
Fredette,  and  D.  N.  Beal  (2008),  “Synchronization  of  Animal-Inspired  Multiple  Fins  in  an 
Underwater  Vehicle  Using  Olivo-Cerebellar  Dynamics,”  IEEE  Journal  of  Oceanic- 
Engineering.  vol.  33,  no.  4,  pp.  563-578. 

Blake.  J.  R.  (1971)  “A  Spherical  Envelope  Approach  to  Ciliary  Propulsion.”  Journal  of  Fluid 
Mechanics,  vol.  46,  part  1.  pp.  199-208. 

Bouzarth,  E.  L.,  A.  Brooks.  R.  Camassa,  H.  Jing,  T.  J.  Leiterman,  R.  M.  McLaughlin,  R. 

Superfine,  J.  Toledo,  and  L.  Vicci  (2007),  “Epicyclic  Orbits  in  a  Viscous  Fluid  About  a 
Precessing  Rod:  Theory  and  Experiments  at  the  Micro-  and  Macro-Scales.”  Physical 
Review  E  Statistical,  Nonlinear,  and  Soft  Matter  Physics,  vol.  76  (part  1  of  2),  016313. 


39 


Brokaw,  C.  J.  (1966),  “Bend  Propagation  Along  Flagella,”  Nature  (London),  vol.  209,  p.  161. 

Burger,  C.,  B.  S.  Hsiao,  and  B.  Chu  (2006),  “Nanofibrous  Materials  and  Their  Applications,” 
Annual  Review  of  Materials  Research,  vol.  36.  pp.  333-368. 

Cordova,  N.  J.,  B.  Ermentrout.  and  G.  F.  Oster  (1992),  “Dynamics  of  Single-Motor  Molecules: 
The  Thermal  Ratchet  Model,”  Proceedings  of  the  National  Academy  of  Sciences, 
Biophysics,  vol.  89,  pp.  339-343. 

Dean.  P..  J.  Porril.  C-F.  Ekerot,  and  H.  Jomteli  (2010),  “The  Cerebellar  Microcircuit  as  an 
Adaptive  Filter:  Experimental  and  Computational  Evidence,”  Nature  Reviews: 
Neuroscience,  vol.  1 1.  pp.  30-43. 

Deitmer,  J.  W.,  H.  Machemer.  and  B.  Martinac  (1984),  “Motor  Control  in  Three  Types  of  Ciliary' 
Organelles  in  the  Ciliate  Stylonychia,”  Journal  of  Comparative  Physiology  A,  vol.  154,  pp. 
113-120. 

dePeyer,  J.,  and  H.  Machemer  (1978),  “Are  Receptor-Activated  Ciliary'  Motor  Responses 
Mediated  Through  Voltage  or  Current?”  Nature,  vol.  276,  pp.  285-287. 

dePeyer.  J.  E..  and  H.  Machemer  (1983),  "Threshold  Activation  and  Dynamic  Response  Range 
of  Cilia  Following  Low  Rates  of  Membrane  Polarization  Under  Voltage-Clamp,”  Journal 
of  Comparative  Physiology’,  vol.  1 50.  pp.  223-232. 

Dillon.  R.  H.,  and  L.  J.  Fauci  (2000),  “An  Integrative  Model  of  Internal  Axoneme  Mechanics  and 
External  Fluid  Dynamics  in  Ciliary  Beating.”  Journal  of  Theoretical  Biology,  vol.  207.  pp. 
415-430. 

Dute.  R..  and  C.  Kung  (1978),  “LUtrastructure  of  the  Proximal  Region  of  Somatic  Cilia  in 
Paramecium  Tetraurelia,”  Journal  of  Cell  Biology,  vol.  78.  pp.  451-464. 

Ellington,  C.  P.  (1999),  “The  Novel  Aerodynamics  of  Insect  Flight:  Applications  to  Micro-Air 
Vehicles,”  Journal  of  Experimental  Biology,  vol.  202,  pp.  3439-3448. 

Griffith.  A.  A.  (1920),  “The  Phenomena  of  Rupture  and  Flow  in  Solids,”  Philosophical 
Transactions  of  the  Royal  Society,  vol.  A221.  pp.  163-198. 

Gueron,  S..  and  K.  Levit-Gurevich  (2001),  “A  Three-Dimensional  Model  for  Ciliary  Motion 
Based  on  the  Internal  9  +  2  Structure,”  Proceedings  of  the  Royal  Society  of  London  B,  vol. 
268,  pp.  599-607. 

Gueron.  S..  and  K.  Levit-Guerevich  (1999a),  “Energetic  Considerations  of  Ciliary  Beating  and 
the  Advantage  of  Metachronal  Coordination,”  Proceedings  of  the  National  Academy  of 
Sciences,  vol.  96,  no.  22,  pp.  12240-12245. 


40 


Gueron,  S..  and  K.  Levit-Guerevich  (1999b),  “A  Three-Dimensional  Model  of  Ciliary  Motion 
Based  on  Internal  9+2  Structure/’  Proceedings  of  the  Royal  Society .  vol.  96,  no.  22,  pp. 
12240-12245. 

Gueron,  S.,  K.  Levit-Gurevich.  N.  Liron.  and  J.  J.  Blum  (1997),  “Cilia  Internal  Mechanism  and 
Metachronal  Coordination  as  the  Result  of  Hydrodynamical  Coupling,"  Proceedings  of  the 
National  Academy  of  Sciences,  Applied  Mathematics ,  vol.  94,  pp.  6001-6006. 

Gueron,  S..  and  N.  Liron  (1992),  “Ciliary  Motion  Modeling,  and  Dynamic  Multicilia 
Interactions,”  Biophysics  Journal .  vol.  63,  pp.  1045-1058. 

Guevorkian.  K..  and  J.  M.  Valles  (2006).  “Swimming  Paramecium  in  Magnetically  Simulated 
Enhanced,  Reduced,  and  Inverted  Gravity  Environments,”  Proceedings  of  the  National 
Academy  of  Sciences,  vol.  103.  no.  35.  pp.  13051-13056. 

Guirao,  B.  (and  fifteen  others)  (2010),  “Coupling  Between  Hydrodynamic  Forces  and  Planar 
Cell  Polarity  Orients  Mammalian  Motile  Cilia."  Nature  Cell  Biology ,  vol.  12,  pp.  341-350. 

Gittes.  F..  B.  Mickey.  J.  Nettleton.  and  J.  Howard  (1993).  “Flexural  Rigidity  of  Microtubules  and 
Actin  Filaments  from  Thermal  Fluctuations  in  Shape,”  Journal  of  Cell  Biology ,  vol.  120, 
no.  4.  pp.  923-934. 


Hines.  M.,  and  J.  J.  Blum  (1978),  “Bend  Propagation  in  Flagella.  1.  Derivation  of  Equations  of 
Motion  and  Their  Simulation,”  Biophysics  Journal,  vol.  23.  pp.  41-57. 

Hines.  M.,  and  J.  J.  Blum  (1983),  “Three-Dimensional  Mechanics  of  Eukary  otic  Flagella," 
Biophysics  Journal,  vol.  41,  pp.  67-79. 

Hohman.  M.  M.,  M.  Shin,  G.  Rutledge,  and  M.  P.  Brenner  (2001),  “Electrospinning  and 

Electrically  Forced  Jets,  II.  Applications,”  Physics  of  Fluids,  vol.  13,  no.  8,  pp.  2221-2236. 

Irwin,  G.  R.  (1948),  “Fracture  Dynamics,”  in  Fracturing  of  Metals,  American  Society  for  Metals, 
Cleveland,  OH.  pp.  147-166. 

Jaffe,  L.  F.  (2007),  “Models  and  Speculations  Stretch-Activated  Calcium  Channels  Relay  Fast 
Calcium  Waves  Propagated  by  Calcium-Induced  Calcium  Influx,”  Biological  Cell ,  vol.  99, 
pp.  175-184  (printed  in  Great  Britain.  doi:10.1042/BC20060031  Scientiae  forum). 

Karp  G.  (2005),  Cell  and  Molecular  Biology:  Concepts  and  Experiments,  Fourth  Edition,  John 
Wiley  and  Sons,  Hoboken,  NJ,  pp.  346-358. 

Kazantsev,  V.  B.,  V.  Nekorkin,  V.  Makarenko,  and  R.  Llinas  (2003),  “Olivo-Cerebellar  Cluster- 
Based  Universal  Control  System,”  Proceedings  of  the  National  Academy  of  Sciences,  vol. 
100.  no.  22,  pp.  13064-13068. 


41 


Kazantsev,  V.  B.,  V.  Nekorkin.  V.  Makarenko,  and  R.  Llinas  (2004),  “Self-Referential  Phase 
Reset  Based  on  Inferior  Olive  Oscillator  Dynamics.  Proceedings  of  the  National  Academy 
of  Sciences,  vol.  101.  pp.  18183-18188. 

Khalil.  H.  K.  (1996),  Nonlinear  Systems,  Prentice-Hall,  Upper-Saddle  River.  NJ. 

Kongthon.  J.,  B.  McKay,  D.  lamratanakul,  K.  Oh,  J.-H.  Chung.  J.  Riley,  and  S.  Devasia  (2010), 
“Added-Mass  Effect  in  Modeling  of  Cilia-Based  Devices  for  Microfluidic  Systems.”  ASME 
Journal  of  Vibration  and  Acoustics,  vol.  132,  024501-1. 

Lauder,  G.V.  (2009).  Private  Communication  with  P.  Bandyopadhyay.  Naval  Undersea  Warfare 
Center  Division,  Newport.  Rl. 

Lauder,  G.  V..  and  E.  G.  Drucker  (2004),  “Morphology  and  Experimental  Hydrodynamics  of 
Fish  Fin  Control  Surfaces,”  IEEE  Journal  of  Oceanic  Engineering ,  vol.  29,  pp.  556-571. 

Lauder.  G.  V.,  P.  G.  A.  Madden.  J.  L.  Tangora.  E.  Anderson,  and  T.  V.  Baker  (2011), 

“Bioinspiration  from  Fish  for  Smart  Material  Design  and  Function,”  Smart  Materials  and 
Structures ,  vol.  20.  094014. 

Li.  F..  and  X.  Yin  (2009),  “Axisymmetric  and  Non-Axisymmetric  Instability  of  an  Electrified 
Viscous  Coaxial  Jet,”  Journal  of  Fluid  Mechanics ,  vol.  632.  pp.  199-225. 

Llinas.  R.  R.  (2009).  “Inferior  Olive  Oscillation  as  the  Temporal  Basis  for  Motricity  and 

Oscillatory  Reset  as  the  Basis  for  Motor  Error  Correction,”  Neuroscience,  vol.  162,  no.  3, 
pp.  797-804. 

Llinas,  R.  R.  (2001),  I  of  the  Vortex  from  Neurons  to  Self  MIT  Press,  Cambridge.  MA. 

Llinas.  R.  R.,  E.  Leznik,  and  V.  I.  Makarenko  (2004),  "The  Olivo-Cerebellar  Circuit  as  a 

Universal  Control  System,”  IEEE  Journal  of  Oceanic  Engineering,  vol.  29,  pp.  631-639. 

Llinas.  R.  R.,  and  Y.  Yarom  (1981),  "Electrophysiology  of  Mammalian  Inferior-Olive  Neurons 
in  Vitro.  Different  Types  of  Voltage  Dependences  in  Conductances,”  Journal  of  Physics 
(London),  vol.  315,  pp.  549-567. 

Lu.  X.  F.,  C.  Wang,  and  Y.  Wei  (2009),  “One-Dimensional  Composite  Nanomaterials:  Synthesis 
by  Electrospinning  and  Their  Applications,  Small,  vol.  5,  pp.  2349-2370. 

Machemer.  H.  (1972),  “Ciliary’  Activity  and  the  Origin  of  Metachrony  in  Paramecium:  Effects  of 
Increased  Viscosity,”  Journal  of  Experimental  Biology ,  vol.  57,  pp.  239-259. 

Machemer,  H..  and  A.  Ogura  (1979),  “Ionic  Conductances  of  Membranes  in  Ciliated  and 
Deciliated  Paramecium.”  Journal  of  Physiology',  vol.  296,  pp.  49-60. 


42 


Machemer.  H.,  and  S.  Machemer-Rehnisch  (1984),  “Mechanical  and  Electric  Correlates  of 

Mechanoreceptor  Activation  of  the  Ciliated  Tail  in  Paramecium.”  Journal  of  Sensory  and 
Comparative  Physiology  A.  vol.  154.  pp.  273-278. 

McCullen.  S.  D.,  K.  L.  Stano,  I).  R.  Stevens.  W.  A.  Roberts.  N.  Monteiro-Riviere,  L.  1.  Clarke, 
and  R.  E.  Gorga  (2007),  “Development,  Optimization,  and  Characterization  of  Electrospun 
Poly(Lactic  Acid)  Nanofibers  Containing  Multi-Walled  Carbon  Nanotubes,”  Journal  of 
Applied  Polymeric  Science,  vol.  105.  pp.  1668-1678. 

Nakaoka.  Y„  H.  Tanaka,  and  F.  Oosawa  ( 1 984),  “Ca~  -Dependent  Regulation  of  Beat  Frequency 
of  Cilia  in  Paramecium,”  Journal  of  Cell  Science ,  vol.  65.  pp.  223-231 . 

Ojha.  S.S.,  D.  R.  Stevens.  K.  Stano,  T.  Hoffman,  W.  E.  Krause,  L.  I.  Clarke,  and  R.  E.  Gorga 
(2008),  “Characterization  of  Electrical  and  Mechanical  Properties  for  Coaxial  Nanofibers 
with  Poly-Ethylene  Oxide  (PEO)  Core  and  Multiwalled  Carbon  Nanotube/PEO  Sheath.” 
Macromolecules ,  vol.  41.  pp.  2509-2513. 

Omoto.  C.  K..  and  C.  Rung  (1980),  “Rotation  and  Twist  of  the  Central-Pair  Microtubules  in  the 
Cilia  of  Paramecium.”.  Journal  of  Cell  Biology,  vol.  87,  pp.  33-46. 

Pemberg.  J„  and  H.  Machemer  (1995),  “Voltage-Dependence  of  Ciliary  Activity  in  the  Ciliate 
Didinium  Nasutum.”  Journal  of  Experimental  Biolog}',  vol.  198,  pp.  2537-2545. 

Rallison.  J.  M.,  and  E.  J.  Hinch  (1995),  “Instability  of  a  High-Speed  Submerged  Elastic  Jet,” 
Journal  of  Fluid  Mechanics ,  vol.  288.  pp.  3 1 1-324. 

Ramaswamy,  S.,  L.  I.  Clarke,  and  R.  E.  Gorga  (201 1),  “Morphological,  Mechanical,  and 

Electrical  Properties  as  a  Function  of  Thermal  Bonding  in  Electrospun  Nanocomposites,” 
Polymer ,  vol.  52.  pp.  3 1 83-3 1 89. 

Ramia,  M..  I).  L.  Tullock.  and  N.  Phan-Thien  (1993),  “The  Role  of  Hydrody  namic  Interaction  in 
the  Locomotion  of  Microorganisms,”  Biophysical  Journal,  vol.  65.  pp.  755-778. 

Reneker.  I).  H..  and  I.  Chun  (1996),  “Nanometre  Diameter  Fibres  of  Polymer  Produced  by 
Electrospinning,”  Nanotechnology,  vol.  7.  pp.  216-223. 

Reneker,  D.  H.,  A.  L.  Yarin.  E.  Zussman.  and  H.  Xu  (2007),  “Electrospinning  of  Nanofibers 
from  Polymer  Solutions  and  Melts,”  Advances  in  Applied  Mechanics,  vol.  41,  pp.  43-195. 

Shin,  Y.  M.,  M.  M.  Hohman,  M.  P.  Brenner,  and  G.  C.  Rutledge  (2001),  “Electrospinning:  A 
Whipping  Fluid  Jet  Generates  Submicron  Polymer  Fibers,”  Applied  Physics  Letters,  vol. 
78.  pp.  1149-1151. 

Shields,  A.  R.,  B.  L.  Fiser,  B.  A.  Evans,  M.  R.  Falvo,  S.  Washburn,  and  R.  Superfine  (2010), 
“Biomimetic  Cilia  Arrays  Generate  Simultaneous  Pumping  and  Mixing  Regimes,” 
Proceedings  of  the  National  Academy  of  Sciences,  vol.  107.  no.  36,  pp.  15670-15675. 


43 


Slotine,  J.-J.  E..  and  W.  Li  (1991),  Applied  Nonlinear  Control ,  Prentice-Hall,  Englewood  Cliffs, 
NJ,  pp.  160-163. 

Tamm.  S.  L.  (1984),  “Mechanical  Synchronization  of  Ciliary  Beating  within  Comb  Plates  of 
Tenophores,”  Journal  of  Experimental  Biolog}’,  vol.  1 13,  pp.  401-408. 

Tamm.  S.  L..  and  S.  Tamm  (1989),  “Calcium  Sensitivity  Extends  the  Length  of  ATP-Reactivated 
Ciliary  Axonemes,"  Proceedings  of  the  National  Academy  of  Sciences,  Cell  Biology ,  vol. 
86,  pp.  6987-6991. 

Teunis,  P.  M.  F.,  and  H.  Machemer  (1994),  “Analysis  of  Three-Dimensional  Ciliary  Beating  by 
Means  of  High-Speed  Stereomicroscopy,’'  Biophysical  Journal,  vol.  67,  pp.  381-394. 

Triantafyllou,  M.  S.,  and  G.  S.  Triantafyllou  (1995),  “An  Efficient  Swimming  Machine,” 
Scientific  American,  vol.  272.  no.  3.  pp.  64-70. 

Yeh.  Y.,  R.  J.  Baskin,  S.  Shen.  and  M.  Jones  (1990).  “Photon  Correlation  Spectroscopy  of  the 
Polarization  Signal  of  Single  Muscle  Fibers,"  Journal  of  Muscle  Research  and  Cell 
Motility,  vol.  1 1,  no.  2.  pp.  137-146. 

Wakeling,  J.  M.,  and  C.  P.  Ellington  (1997),  “Dragonfly  Flight.  III.  Lift  and  Power 
Requirements,”  Journal  of  Experimental  Biolog}’,  vol.  200,  pp.  583-600. 

Wurzel,  S.  E.  (2003),  “Prospects  for  Studying  Negative  Gravitaxis  in  Paramecium  Using 
Magnetic  Field  Gradient  Levitation,”  Senior  Thesis.  Brown  University.  Department  of 
Physics,  Providence,  RI. 

Zwillinger,  I).  (1996),  Standard  Mathematical  Tables  and  Formulae,  CRC  Press. 


44 


APPENDIX— SUPPLEMENTARY  INFORMATION 


A.l  TWO-DIMENSIONAL  HYDRODYNAMIC  MODELING 


Note:  This  section  is  due  to  Norman  Toplosky  (Code  1532). 


As  a  first  step  in  modeling  the  dynamics  of  cilium  motion,  a  two-dimensional  dynamic 
model  was  constructed  to  develop  familiarity  with  the  solution  technique.  Following  the 
biorobotic  approach  mentioned  earlier,  the  authors  made  no  direct  effort  at  modeling  the  actual 
physiology  of  the  cilium.  but  instead  concentrated  on  driving  the  base  of  a  mechanically- 
equivalent  cilium  with  the  same  motion  as  observed  in  the  data  of  Teunis  and  Machemer  (1994), 
and  then  attempted  to  adjust  the  rigidity  and  drag  parameters  so  as  to  reproduce  the  motion  of  the 
entire  cilium.  The  two-dimensional  model  was  considered  to  be  an  important  first  step  and 
building  block  toward  the  eventual  development  of  a  full  three-dimensional  model.  The  latter 
model  would  be  different  from  the  former  in  two  important  ways;  first,  a  third  dimension  would 
allow  out-of-plane  bending  and  torsion,  and  second,  the  three-dimensional  model  would  possibly 
allow  inclusion  of  torsional  rigidity  and  a  driving  torque  at  the  cilium  base. 

In  the  initial  two-dimensional  model,  the  cilium  is  modeled  as  A  rigid  links  connected  by 
nodal  ball  joints  in  which  the  cilium  bending  rigidity  is  lumped.  The  bending  rigidity  lumped  at 
the  nodes  can  vary  both  with  time  (depending  on  where  one  is  in  the  power-recovery  stroke 
cycle)  and  with  length  along  the  cilium.  Fluid  drag  loads — due  both  to  the  mean  advance  of  the 
cilium  (velocities  U,  V)  and  to  the  rotational  motion  of  the  cilium  about  its  base — act  on  each 
link  and  depend  on  both  the  orientation  of  the  link  (*.>’,  and  z  coordinates  of  the  bracketing 
nodes)  and  the  in-water  velocity  of  the  link  (u,  v.  and  w  velocity  components  of  the  bracketing 
nodes).  The  bending  moment  in  the  nodal  ball  joints  is  proportional  to  the  local  curvature, 
estimated  by  the  change  in  orientation  of  the  tw  o  links  joined  at  that  node.  Internal  forces  also 
act  at  each  node.  Figure  A-l  shows  this  setup.  Note  that  in  the  biological  cilium  and  the  robotic 
cilium.  y  represents  the  vertical  direction  and  z  represents  the  spanwise  direction;  but  in  the  two- 
dimensional  hydrodynamic  simulation,  z  represents  vertical  direction. 

The  polar  angle  of  the  first  link  (the  angle  between  the  link  and  the  vertical  r-axis)  and  its 
rate  (#and  dOldt)  are  taken  as  kinematic  boundary  conditions  at  the  base.  Zero  force  and 
moment  conditions  are  applied  at  the  distal  node  of  the  last  link. 

The  physics-based  two-dimensional  model  using  N  links  is  well-posed:  Not  counting  the 
origin  (cilium  base  is  node  0),  there  are  Anodes  and  at  each  of  these  there  are  two  position 
coordinates,  two  velocity  components,  and  two  force  components,  for  a  preliminary  total  of  6N 
unknowns  (see  table  A-l).  On  each  of  the  N  links  there  are  two  force  balances  (in  the  x-  and  z- 
directions).  one  moment  balance  (in  the  y-direction).  one  geometric  constraint  ensuring  that  the 
link's  bracketing  nodal  positions  are  a  link  length  apart,  and  two  kinematic  relations  relating  the 
distal  node’s  velocity  components  to  the  time  derivatives  of  the  position  coordinates,  for  a 
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preliminary  total  of  6N  equations.  However,  this  preliminary  set  of  6N  equations  in  6 N 
unknowns  is  slightly  modified  by  the  additional  constraints  and  unknowns  at  the  base  and  distal 
node. 

The  resulting  equations  (see  table  A-l)  were  programmed  in  MATLAB.  The  base  link  was 
driven  by  a  harmonic  motion  (the  polar  angle  was  sinusoidal).  Cilium  stiffness  was  increased 
during  the  power  stroke  and  relaxed  during  the  recovery'  stroke,  and  a  mean  speed-of-advance  (as 
a  fraction  of  the  cilium  tip  speed  midway  through  the  power  stroke)  was  added  to  the  base 
motion.  The  profiles  in  figure  A-2  were  computed. 

The  development  of  the  three-dimensional  cilium  model  is  still  in  progress.  As  mentioned 
earlier,  there  are  two  major  changes  when  going  from  two  dimensions  to  three.  The  first  change 
is  that  all  vectors  have  three  components  instead  of  two.  and  because  bending  is  no  longer 
confined  to  the  y-direction,  the  cilium  develops  torsion  as  the  bending  direction  (the  unit 

bi-normal  vector  b  )  changes  direction.  The  second  change  is  that  one  can  include  a  driving  base 
torque  and  resulting  twist,  if  one  wishes. 

Table  A-l  provides  background  information  regarding  the  two-dimensional  hydrodynamic 
modeling.  This  section  also  lists  the  Frenet  equations,  which  describe  the  geometry'  of  a  curve  in 

space;  here,  k  is  curvature,  r  is  torsion,  t  and  fi  are  unit  tangent  and  normal  vectors,  b  is  the 
unit  bi-normal  vector,  and  s  is  the  axial  length. 
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Figure  A-2.  Two-Dimensional  Ciliunt  Profiles  (nu:  normalized  unit) 

The  integral  of  the  base  force  over  a  beat  cycle  is  the  net  impulse  from  that  cycle.  The 
asymmetric  beating  pattern  and  mean  advance  of  the  base  allows  for  a  non-zero  thrust,  even  in 
two  dimensions,  as  shown  in  figure  A-3. 
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Figure  A-3.  Two-Dimensional  Base  Thrust  (nu:  normalized  unit) 
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Table  A-l.  Hydrodynamic  Modeling  Background  Information 


Unknowns 

Equations 

Node  0  (base) 

Fx 

Fz 

My 

Node  1 

X 

x  =  Ls\n0 

z 

z  =  LcosO 

u 

u  =  L(dO/dt)  cos  6 

V 

Fx 

v  =  -L(d0/dt)s\nO 

Fz 

Nodes  2... N-  1 

X 

U=  dx/dt 

z 

V  =  dy/dt 

u 

V 

Fx 

Fz 

Node  N  (distal  node) 

X 

U=  dx/dt 

z 

V=  dy/dt 

u 

V 

Link  1 

Fx  force  balance 

Fy  force  balance 

My  moment  balance 

Link  2 ...N 

Fx  force  balance 

Fy  force  balance 

My  moment  balance 
Link-length  constraint 

Total: 

67V+1 

6  AM- 1 

Frenet  Equations: 


di 

—  =  +Kn 
ds 


d  n 
ds 


=  -Kt  +  Tb 


dh 

—  =  -r  n 
ds 
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A.2  VALIDATION  OF  DIGITIZATION  OF  PUBLISHED  BIOLOGICAL 
CILIUM  DATA 

A.2. 1  Parametric  Fits 

The  digitized  biological  cilium  in  pixels  in  the  longitudinal  plane  is  shown  in  figure  A-4. 

The  steps  for  obtaining  the  three-dimensional  parametric  fits  and  associated  lengths  are  as 
follows: 

1 .  Find  the  base  of  the  top  and  side  views  and  adjust  the  pixel  values  for  a  base  of  (0,  0). 

2.  Obtain  the  length  in  the  orthogonal  axis  (from  the  alternate  view)  for  the  top  and  side 
projections  (they-value  in  the  top  view  is  z  (width)  and  the  y- value  in  the  side  view'  is  >•  (height)): 

a.  for  the  side  view,  this  is  the  top  view'  vertical  axis. 

b.  for  the  top  view,  this  is  the  side  view  vertical  axis. 

3.  Interpolate  the  orthogonal  lengths  (from  above)  to  match  the  number  of  points  from  the 
alternate  view.  This  is  required  since  the  top  view  and  side  view  of  a  particular  time  step  have  a 
different  number  of  points. 

4.  Calculate  the  lengths  for  each  view  by  V(Ax:+  Ay:+  Az2).  where  x  and  y  are  from  the  top 
and  side  planar  views  and  z  is  the  orthogonal  axis  differences  from  step  3.  Note  that  inconsistent 
data  between  the  two  views  result  in  the  lengths  for  each  view  being  different. 

5.  Obtain  two-dimensional  parametric  fits  for  the  top  and  side  views  using  the  lengths 
calculated  above,  with  the  length  as  the  independent  variable. 

6.  Choose  the  x ,  y,  z  values  (to  use  in  the  three-dimensional  parametric  fit  later)  as  follows: 

a.  x  =  top  view  x  (could  also  choose  side  view  x  values). 

b.  y  =  side  view  y. 

c.  z  =  top  viewy. 

7.  Calculate  the  lengths  as  the  cumulative  sum  of  V(Ax2+  Ay2+  Az2)  for  each  time  step.  See 
figure  3a  in  the  main  text. 

8.  Use  the  parametric  equations  for  both  views  to  calculate  an  arbitrary,  user-selected 
number  of  orbits  for  analysis. 

9.  Fit  a  three-dimensional  parametric  equation  to  the  data  for  each  orbit  using  the  .y,  y,  z 
values,  with  the  step  number  as  the  independent  variable. 

10.  Calculate  the  lengths  as  the  cumulative  sum  of  V(Ax2+  Ay2+  Az2)  for  each  time  step, 
through  each  orbit.  See  figure  3  c. 
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Figure  A-4.  Graph  of  Digitized  Citium  Prior  to  Parametric  Processing  and  Smoothing 
(Each  color  is  a  different  time  step;  length  units  are  in  pixels.) 


A.2.2  Notes  on  Data  Smoothing  and  Definitions 

1 .  A  length  of  1 56  is  the  original  pixel  data  (digitized  biological  data);  a  length  of  1 7  pm  is 
for  the  adjusted  biological  cilium;  and  a  length  of  0.21  m  is  for  the  robotic  cilium. 

2.  The  polynomial  used  along  the  length  of  the  cilium  is  fourth  order,  although  it  is  used 
only  to  calculate  an  equal  number  of  points  for  each  time  step. 

3.  The  polynomial  for  each  orbit  (concentric  path)  is  15th  order. 

4.  Parametric  equations  are  independent  for  x,  y,  z  and  are  based  only  on  length  and 
themselves,  i.e.,  x  =fx.l).y  =J(y,l),  z  =J{z.l),  where  /  is  the  length  of  the  cilium  (/  =  5). 
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5.  Smoothing  is  done  with  a  five-point  moving  average  to  obtain  the  average  cilium  radius 
(curved  length  projected  on  the  transverse  plane y-z)  at  a  given  time  instant: 

w  here  /  (/)  is  the  curved  length  of  the  entire  cilium  in  the>’-r  plane  (  transverse  plane)  (sum  of 

the_y-  and  --components  of  each  length  segment  of  the  cilium)  at  a  given  time  instant;  to  clarify, 
Ravg  is  not  a  radial  distance  from  the  origin  to  the  point  on  the  cilium  where  Ring  lies. 

6.  The  accuracy  of  the  Teunis-Machemer  data  digitization  and  length  calculation  is  as  given 
in  tables  A-2  through  A-4.  The  cilium  length  is  divided  into  20  segments  and  the  lengths  are 
measured  from  the  base. 


A. 2.3  Accuracy  of  Parametric  Fits  and  Smoothing 

The  accuracy  of  the  parametric  fits  and  smoothing  is  shown  in  tables  A-2  through  A-4. 
These  tables  show  how  the  mean  and  mis  distances  (pm)  of  20  points  on  the  cilium  counting 
vary’.  Table  A-2  is  for  the  entire  beat  cycle;  tables  A-3  and  A-4  are  for  the  power  and  return 
strokes,  respectively.  (In  the  column  headings  of  these  tables,  “top”,  “mid”,  and  “hot”  refer  to 
the  top  (a),  middle  (b),  and  bottom  (c)  frames  in  figure  3  in  the  main  text.)  Although  the  cilium 
is  curved  during  the  return  stroke,  the  cilium  lengths  are  close  (within  4%)  to  those  during  the 
power  stoke  when  the  cilium  is  fairly  straight.  The  data  processing,  therefore,  is  accurate. 
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Table  A-2.  Variation  of  RMS  in  Length  Estimation  in  the  Revisited  TM  Data  (Entire  Beat  Cycle) 


Mean  and  RMS  Cilium  Lengths  (All  Positions)  (fun) 

Mean  (top) 

Mean  (mid) 

Mean  (bot) 

RMS  (top) 

RMS  (mid) 

RMS  (bot) 

0.1128 

0.0875 

0 

0.1271 

0.0940 

0 

1.0402 

1.0026 

0.9157 

1.0465 

1.0073 

0.9193 

1.9508 

1.9021 

1.8157 

1.9573 

1.9071 

1.8199 

2.8514 

2.7915 

2.7053 

2.8572 

2.7958 

2.7091 

3.7456 

3.6739 

3.5878 

3.7504 

3.6774 

3.5909 

4.6355 

4.5513 

4.4651 

4.6395 

4.5541 

4.4676 

5.5220 

5.4249 

5.3386 

5.5257 

5.4273 

5.3408 

6.4064 

6.2962 

6.2097 

6.4098 

6.2982 

6.2116 

7.2900 

7.1665 

7.0799 

7.2932 

7.1685 

7.0818 

8.1746 

8.0379 

7.9511 

8.1777 

8.0397 

7.9529 

9.0621 

8.9121 

8.8253 

9.0653 

8.9139 

8.8271 

9.9546 

9.7910 

9.7041 

9.9578 

9.7928 

9.7060 

10.8531 

10.6755 

10.5887 

10.8565 

10.6774 

10.5908 

11.7578 

11.5658 

1 1 .4792 

11.7616 

11.5680 

11.4815 

12.6675 

12.4608 

12.3744 

12.6718 

12.4633 

12.3772 

13.5795 

13.3578 

13.2717 

13.5846 

13.3609 

13.2751 

14.4896 

14.2529 

14.1672 

14.4955 

14.2566 

14.1712 

15.3925 

15.1409 

15.0555 

15.3991 

15.1452 

15.0602 

16.2822 

16.0154 

15.9305 

16.2893 

16.0202 

15.9357 

17.1533 

16.8702 

16.7857 

17.1605 

16.8752 

16.7911 
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Table  A-3.  Variation  of  RMS  in  Length  Estimation  in  the  Revisited  TM  Data  (Power  Stroke) 


Mean  and  RMS  Cilium  Lengths  (Power  Stroke:  TS  38-14)  (pm) 

Mean  (top) 

Mean  (mid) 

Mean  (hot) 

RMS  (top) 

RMS  (mid) 

RMS  (bot) 

0.1253 

0.1062 

0 

0.1420 

0.1106 

0 

1.0862 

1.0541 

0.9496 

1.0915 

1.0572 

0.9529 

2.0113 

1.9682 

1.8649 

2.0168 

1.9723 

1.8691 

2.9130 

2.8590 

2.7564 

2.9180 

2.8630 

2.7605 

3.8002 

3.7344 

3.6320 

3.8045 

3.7380 

3.6357 

4.6788 

4.6000 

4.4975 

4.6826 

4.6031 

4.5008 

5.5525 

5.4595 

5.3568 

5.5559 

5.4622 

5.3597 

6.4235 

6.3154 

6.2124 

6.4267 

6.3178 

6.2150 

7.2929 

7.1692 

7.0660 

7.2961 

7.1715 

7.0684 

8.1615 

8.0222 

7.9187 

8.1647 

8.0243 

7.9210 

9.0299 

8.8749 

8.7714 

9.0330 

8.8769 

8.7736 

9.8983 

9.7279 

9.6246 

9.9015 

9.7298 

9.6267 

10.7672 

10.5816 

10.4787 

10.7703 

10.5834 

10.4807 

11.6366 

11.4361 

1  1.3337 

11.6397 

11.4377 

11.3356 

12.5064 

12.2909 

12.1893 

12.5094 

12.2925 

12.1911 

13.3759 

13.1454 

13.0447 

13.3789 

13.1469 

13.0464 

14.2442 

13.9984 

13.8987 

14.2472 

13.9998 

13.9003 

15.1099 

14.8480 

14.7494 

15.1127 

14.8493 

14.7509 

15.9711 

15.6918 

15.5943 

15.9739 

15.6931 

15.5957 

16.8264 

16.5271 

16.4306 

16.8292 

16.5283 

16.4320 
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Table  A-4.  Variation  of  RMS  in  Length  Estimation  in  the  Revisited  TM  Data  (Return  Stroke) 


Mean  and  RMS  Cilium  Lengths  (Return  Stroke:  TS  15-37)  (pm) 

Mean  (top) 

Mean  (mid) 

Mean  (bot) 

RMS  (top) 

RMS  (mid) 

RMS  (bot) 

0.1025 

0.0721 

0 

0.1133 

0.0777 

0 

1.0021  1 

0.9600 

0.8877 

1.0078 

0.9640 

0.8905 

1.9008 

1.8475 

1.7750 

1.9068 

1.8515 

1.7782 

2.8005 

2.7357 

2.6631 

2.8059 

2.7391 

2.6658 

3.7006 

3.6239 

3.5513 

3.7052 

3.6266 

3.5534 

4.5997 

4.5110 

4.4383 

4.6037 

4.5132 

4.4400 

5.4969 

5.3964 

5.3236 

5.5005 

5.3982 

5.3251 

6.3924 

6.2803 

6.2075 

6.3958 

6.2820 

6.2089 

7.2876 

7.1643 

7.0914 

7.2909 

7.1659 

7.0928 

8.1853 

8.0509 

7.9779 

8.1885 

8.0525 

7.9792 

9.0888 

8.9429 

8.8698 

9.0919 

8.9444 

8.8711 

10.0010 

9.8430 

9.7698 

10.0041 

9.8445 

9.7711 

10.9240 

10.7530 

10.6796 

10.9272 

10.7544 

10.6809 

11.8579 

11.6730 

11.5994 

11.8613 

11.6744 

11.6008 

12.8006 

12.6011 

12.5273 

12.8045 

12.6027 

12.5288 

13.7477 

13.5333 

13.4592 

13.7521 

13.5351 

13.4610 

14.6923 

14.4632 

14.3890 

14.6974 

14.4653 

14.3911 

15.6260 

15.3828 

15.3084 

15.6317 

15.3853 

15.3110 

16.5392 

16.2828 

16.2083 

16.5453 

16.2856 

16.2112 

17.4234 

17.1536 

17.0791 

17.4294 

17.1565 

17.0821 
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A. 2.4  Definitions 


In  figure  5,  the  angles  off  the  x-axis  are  defined  as: 
position  angle  =  acos(x  length/total  length)  *  1 80/7T, 


and 


velocity  angle  =  acos(x_velocity/total_velocity)  *  1  80/tt. 


A.2.5  Waveforms  Used  for  Robotic  Control 

In  the  early  phase  of  the  work,  the  data  in  figure  A-4  were  compared  w  ith  a  sine  wave, 
applying  a  phase  difference  to  determine  if  the  biological  cilium  motion  could  be  reproduced  by 
sine  waveform  in  the  robotic  cilium.  Comparison  of  the  angles  with  sine  waves  was  made  by 
adding  a  phase  shift  of  90°  or  120°.  and  the  result  is  shown  in  figure  A-5.  A  sine  wave  is  a  good 
approximation  of  the  biological  cilium  position  and  the  120°  phase  difference  has  a  better  fit  to 
the  biological  cilium  than  the  90°  phase  difference  does.  Both  the  Cartesian  and  the  polar 
coordinates  can  be  used  for  robotic  control. 

The  fit  of  the  sine  wave  suggests  that  the  moment  applied  at  the  base  and  the  motion  of  the 
biological  cilium  are  related,  although  the  cilium  is  alternately  taut  and  slack  during  the  power 
and  return  strokes,  respectively. 
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Figure  AS.  Biological  Cilium  Waveform  from  the  Revisited  TM  Data  (Termed  “ real”  in  Figure 
Legend)  at  the  Distal  Point  Compared  with  Sine  Waveforms  in  Cartesian  and  Polar  Coordinates 
(The  sine  waveforms  (termed  “ abstracted ”)  are  also  compared  with  sine  waves  at  phase  lead/lag 
differences  of  120°  (upper  set  of plots)  and  90°  (lower  set  of  plots).  Left  Column:  Cartesian; 
middle  column:  polar  (real:  biological  cilium;  dark  symbols);  right  column:  polar  (abstracted: 

sine  wave;  light  symbols).) 
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A.3  ROBOTIC  CILIUM 


A.3.  1  Early  Work  on  CAD  Simulation  with  D-Cam  Guide  and  Fabrication  of  a 
D-Cattt  Pendulum 

A.3. 1.1  CAD  Simulation.  Computer-aided  design  (CAD)  is  an  accurately  dimensioned,  three- 
dimensional  geometric  relationship  between  objects  that  is  also  precisely  defined  at  all  instants 
of  time  when  relative  motion  between  the  components  exists.  Hydrodynamic  simulations  of 
living  animals  frequently  lack  fluid-structure  interactions  because  the  structural  properties  are 
unknown.  Therefore,  at  an  early  stage,  it  was  thought  to  be  prudent  to  carry  out  a  CAD 
simulation  to  understand  the  IOOMs  that  are  involved  in  the  cilium  motion  (see  figure  1  la  in  the 
main  text). 

It  is  hypothesized  that  a  D-cam  would  force  the  cilium  to  generate  the  observed  motion 
verifiable  in  orthogonal  planes.  Figure  A-6  shows  successive  frames  from  the  CAD  simulation. 
The  cilium  was  given  planar  bending  only;  hence,  it  will  not  produce  torsion  and  the  question 
arises:  What  is  the  origin  of  torsion?  (see  section  A.4  of  this  appendix).  In  view  of  the 
similarities  of  the  robotic  motion  to  the  biological  cilium.  it  is  concluded  that  the  points  at  any 
fixed  distance  from  the  base  circumscribe  the  shape  given  by  the  letter  D  in  planform — see 
figure  1 3a  (left  column).  The  power  stroke  roughly  extends  over  the  straight  part  of  the  D.  and 
the  return  stroke  traces  the  curved  part.  This  exercise  helped  build  the  cilium  actuation  apparatus 
shown  in  figures  A-7  and  A-8. 


A.3. 1.2  Fabrication  of  a  D-Cam  Pendulum.  Section  A.3. 1 . 1  suggests  that  the  cilium  could  be 
modeled  as  a  compound  pendulum.  Figure  A-7  shows  a  schematic  and  a  photograph  of  the  drive 
and  guide  of  such  a  pendulum  (the  cilium  is  not  shown).  The  positive  and  negative  cams  d+  and 
d-  act  as  the  guide  of  cilium  track.  The  pendulum  base  is  located  at  X.  The  cilium  would  be 
attached  to  the  end  of  C. 

In  figure  A-7.  T,  M.  and  B  are  three  stacked  discs;  disc  B  has  a  cutout  of  shape  D  labeled 
“d-”;  there  is  a  D-shaped  cam  (d+)  between  discs  M  and  T;  a  bar  (L-C)  runs  from  d+  to  d-  at  the 
bottom  of  which  the  cilium  hangs;  there  is  a  universal  support  bearing  X  where  this  bar  meets 
disc  M;  there  is  a  paper  clip  L  that  is  used  by  a  lever  above  the  disc  T  (hidden  under  the  hand  in 
the  photograph)  to  turn  the  C-bar  along  the  path  shown  by  the  broken  line  at  the  top  of  the 
schematic  whereby  the  C-bar  runs  while  hugging  the  two  D-cams  along  the  arrows  shown  or  in 
the  reverse  direction. 

Note  the  similarity  between  figures  A-6  and  A-7 — the  CAD  and  the  hardware  renditions. 
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Figure  A-6.  Successive  Frames  (Time  Increases  from  a  to  d)  from 
CAD  Simulation  of  D-Cam-Guided  Three-Link  Robotic  Cilium 
(The  robotic  cilium  is  straight  in  the  flat  portion  of  the  D-cam, 
and  it  droops  while  in  the  round  part  of  the  cilium.) 


Note:  See  the  operation  of  the  D-cam  in  the  CAD  video  file  titled 
“flagellum  dual.mov'’  in  shared  folder  D:\Documents  and 
Settings\promode.bandyopadhyay\My  Documents\Cilium  (also 
available  from  the  lead  author). 
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Figure  A-7.  Pendulum  Apparatus  for  Understanding  Cilium  Motion: 
(a)  Schematic,  (b)  Photograph  of  the  Apparatus 


A.3.2  Early  Version  of  Robotic  Cilium  Apparatus 

Figure  A-8  provides  a  schematic  and  photograph  of  the  early  version  of  the  robotic  cilium. 
and  figure  A-9  shows  how  planar  bending  is  produced  in  the  robotic  cilium  using  a  pair  of 
antagonistic  sliders. 

In  figure  A-8.  the  symbols  are  as  follows:  U:  universal  joint,  R:  cilium  bending  mechanism 
(  see  figure  A-9),  1-1  and  2-2  make  one  pair  of  bucket  handles,  and  3-3  and  4-4  make  another  pair 
(the  bucket  handles  allow  the  base  of  the  cilium  (U)  to  remain  at  the  same  location  while  the 
IOOM  motors  impart  motion  to  the  cilium);  and  Mx,  M„  and  Mt,  are  x,  z,  and  bending  motors. 
Note  that  in  this  apparatus,  the  twist  motor  (My)  is  absent,  the  need  for  which  became  clear  later. 
This  apparatus  worked,  but  it  was  desirable  to  replace  the  bucket  handle  design  to  reduce 
friction. 


Note:  See  the  video  “FlagModel"  (9,855-kb  Windows  Media; 
available  from  the  lead  author)  to  observe  the  operation  of  the 
dual-strip  cilium  model  in  figure  A-9.  and  of  the  dual  bucket 
handle  mechanism  shown  in  figure  A-8. 
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Figure  A-8.  Schematic  and  Photograph  of  Early  Version  of  the  Biorohotic  Cilium  Apparatus 
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Robotic  Paramecium: 


Slider 


L 


^  -■  The  other  strip  is  glued  at  the  end  and  slider  = 

r'  4  4  fatta 


V  t  t  t  t  t  I 

One  carbon  fiber  strip  is  glued  at  each  ring  and  bas 


i  ring  and  base 


End  view  of 
end  ring 


Figure  A-9.  Schematic  of  Mechanism  Showing  How  Planar  Bending  Is  Produced 
in  the  Robotic  Cilium  by  Sliding  the  Top  Strip  Relative  to  the  Bottom  Strip 
(Both  sketches  are  side  views.) 


A.3.3  Experimental  Simulation  of  Cilium  Reynolds  Number 

Using  catheter  tubing  of  1  mm  diameter  for  a  flow  velocity  U  scale  of  1  mm/s,  and  water  at 
5°C,  the  kinematic  viscosity  is  1.5  x  10'f'  m2/s.  giving  a  Reynolds  number  Re  of  0.66  .which  is  < 
1 .0  for  the  same  drag-versus-/te  relationship  to  apply.  Therefore,  using  a  magneto-rheological 
fluid,  it  may  be  possible  to  simulate  the  property  of  variable  hardness  in  a  realistic  Re  number 
and  perhaps  neutral  buoyancy.  An  even  lower  Reynolds  number  is  desirable  and  this  is 
achievable  for  U=  0.1  mm/s.  which  gives  an  Re  =  0.066  (which  is  <  0.1). 
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A.4  “FRICTIONLESS”  CILIUM  ACTUATION  USING  A  HEMISPHERICAL 
MOTOR  DRIVE 

In  a  pendulum  design,  it  is  important  to  lower  friction.  Figure  A- 10  shows  photographs  of  a 
hemispherical  motor,  and  figure  A-l  1  is  a  schematic  of  the  magnets  and  air  gaps.  This  motor 
should  be  an  extremely  low  friction  pendulum,  which  is  an  alternative  to  the  Maxon  motor  drives 
used  in  figures  1 1  and  A-8. 


\ 

\ 

%  •  ^  # 

f*  ""JL  f*  ^  m 

1 
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Electromagnets 

Stator 

(c) 

(d) 

Figure  A-l  0.  Photographs  of  the  Rotor  (b)  and  Stator  (d)  of  a 
Hemispherical  Motor  (a,  c)  Assembly 
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Figure  A-l  I.  Schematic  of  Electromagnets  and  Gaps  in  the  Hemispherical  Motor 


The  hemispherical  motor  has  35  electromagnets.  The  power  supply  to  these  motors  is  to  be 
programmed  to  drive  the  rotor.  I  he  program  should  ensure  that  the  required  IOOMs  are 
imparted  on  the  cilium.  Figure  A-lOa  is  view  of  the  assembly:  A-lOb  shows  the  underside  of  the 
top  plate  of  the  pendulum  drive  and  the  gimbal:  A- 10c  shows  a  side  view  of  the  assembly:  and 
A-lOd  shows  the  35-electromagnet  assembly  exposed:  the  pendulum  ball  in  figure  A-lOb  skirts 
the  hemispherical  surface  resulting  from  the  assembly  shown  in  A-lOd  with  a  small  clearance. 


Figure  A-l  1  is  a  schematic  of  the  electromagnets  and  of  the  air  and  steel  gaps.  This 
schematic  was  used  to  model  the  force  produced  by  each  electromagnet.  The  following 
formulation  is  due  to  J.  Dana  Hrubes  (NUWC  Code  8233).  who  had  used  this  relationship  and 
measurements  with  a  moment  arm  to  estimate  the  viability  of  the  hemispherical  motor  shown  in 
figures  A- 10  and  A-l  1 .  The  magneto-motive  force  F  is  given  by 


(nQ2  Mo  A 

2£  2 

t&P 


(A-l) 
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where  n  =  the  number  of  turns  in  the  coil  in  each  electromagnet,  /  =  the  current  through  the  coil, 
ju0  =  the  magnetic  permeability  in  a  vacuum,  A  =  area,  and  £  =  the  air  gap  between  the 

magnet  and  steel  in  the  direction  of  the  magnetic  flux.  Compared  to  ju0 ,  the  magnetic 
permeability  of  cold  rolled  steel  is  2000  times  larger.  For  a  gap  £  of  3  to  4  mm  at  /  =  1 .5  A, 

both  equation  (A-l)  and  measurement  of  maximum  force  lead  to  a  force  of  0.74  N.  The  tentative 
conclusion  is  that  the  design  in  figure  A- 10  would  produce  the  required  cilium  motion. 


A.5  ROLE  OF  TORSIONAL  RIGIDITY  IN  THE  MANIFESTATION  OF  TORSION 

The  objective  is  to  understand  the  role  of  material  properties  (such  as  torsional  elasticity  and 
moment  of  inertia)  in  the  manifestation  of  torsion  (production  of  off-plane  curv  ature)  in  a  surface 
under  fluid  loading.  Several  three-doublet-long  robotic  cilia  have  been  fabricated.  As  figure 
A- 12  shows,  the  three  doublets  consist  of  three  pairs  of  fixed-length,  inextendable  tension  wires 
positioned  diametrically  across  within  the  pair,  and  120°  apart  between  pairs;  within  a  pair,  the 
tension  wire  is  continuous  and  is  given  two  windings  over  the  pulley.  There  is  also  a  spring- 
loaded  central  tension  wire.  The  tension  wires  pass  through  the  cilium,  which  is  made  of 
segments  of  Plexiglas  rods.  There  is  a  gap  of  about  1  mm  betw  een  these  rod  segments;  a  small 
ball,  which  acts  as  a  pivot,  is  positioned  at  the  center  between  adjacent  rod  segments.  A  real 
cilium  has  9+2  doublet  (“axenome'”)  construction  (Dute  and  Rung,  1978).  Like  the  biological 
doublets,  the  robotic  doublets  also  produce  relative  motion  (sliding)  between  a  pair  of  axial  bars. 
Three-dimensionality  develops  in  the  abstracted  robotic  cilium  in  figure  A- 12  as  follows: 

When  neutrally  buoyant,  the  geometric  deformation  ( k .  r)  in  the  robotic  cilium  is  given  by 

(k,  x)=fn(G,E,I,Cd,T), 

w  here  G  and  E  are  the  moduli  of  torsion  and  elasticity,  respectively,  /  is  the  moment  of  inertia, 
Cd  is  the  coefficient  of  drag,  and  T  is  the  applied  distortion  torque. 

In  the  particular  example  in  figure  A- 13.  when  a  torque  is  applied  to  one  of  the  levers,  a 
helical  spiral  is  produced  (curvature  k  and  torsion  r  are  both  constant)  and  the  deflection  is  not 
planar.  This  finding  came  as  a  surprise,  because  a  planar  distortion  (as  shown  in  figure  A- 12) 
was  expected.  The  levers  were  found  to  be  coupled  (if  one  is  rotated,  the  others  also  rotate).  So, 
the  question  was:  What  is  the  mechanism  of  the  three-dimensional  distortion  produced? 

As  figure  A-l 4  shows,  when  sheathing  is  tightly  slipped  over  the  three-doublet  cilium, 
planar  deflection  is  produced;  that  is.  a  helical  spiral  is  not  produced.  (The  robotic  cilium  in 
figure  A- 14  is  shorter  than  that  shown  in  figure  A-l  3.)  The  torsion  does  not  manifest  itself  this 
time.  This  behavior  can  be  understood  by  recognizing  the  role  of  the  modulus  of  torsional 
rigidity.  Note  that  (as  figure  A- 12  shows)  the  robotic  cilium  has  air  gaps  between  the  segments. 
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Figure  A-I2.  Schematic  of  Mechanism  for  Actuating  a  Long  Robotic  Cilium 


(a)  (b) 


Figure  A-13.  Development  of  Torsion  in  a  Low-Hardness  “Rod” 
(In  (a)  and  (b),  the  lower  thumbwheel  is  being  turned  in  the  clockwise 
direction.  Broken  lines:  Lever  and  baseline  robotic  cilium  axes. 
Deformation  is  propagated  from  the  distal  point  to  the  base.) 
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Figure  A-14.  Effects  of  Clockwise  and  Anti-Clockwise  Turning  of  a  Doublet  of  a  Two-Doublet 
Sheathed  Cilium  Oriented  Orthogonally:  (a)  and  (c)  -  Doublet  1;  (d)  and  (e)  -  Doublet  2 
(Broken  lines:  Lever  and  baseline  robotic  cilium  axes.) 


A-22 


which  makes  the  torsional  rigidity  of  the  cilium  extremely  low  (close  to  zero).  On  the  other 
hand,  when  tight  sheathing  is  slipped  over  it,  the  cilium  becomes  hard.  Therefore,  high  hardness 
resists  the  manifestation  of  torsion ,  hut  low  hardness  does  not.  Importantly,  a  critical  hardness 
should  exist  at  the  boundary  of  the  two  hardnesses  in  the  hardness  versus  torsion  relationship, 
where  a  slight  increase,  or  decrease  in  hardness,  would  produce  substantial  change  in  the 
manifestation  of  torsion. 

Let  if/  be  the  angular  distortion/axial  length.  Then, 

GI  =  7>. 

If  y/  remains  the  same  in  samples  1  and  2, 

TJ{GI\=Tf{GI\. 

For  the  given  angular  deflection,  when  GI  is  low,  less  torque  would  be  required  to  shear  to  the 
same  angle  over  the  same  axial  distance.  (When  G  is  reduced  from  a  high  value,  the  levers  in 
figure  A- 13  are  easier  to  turn.)  Therefore,  GI  should  be  as  low  as  possible  so  that  the  work 
required  to  deform  is  reduced. 

The  sensitivity  of  resulting  torsion  and  applied  torque  to  hardness  w  as  further  verified  in  the 
following  manner.  Two  robotic  rubber  cilia  (figure  A- 15)  of  two  different  hardnesses  were 
fabricated  (the  segment  gaps  shown  in  figure  A- 12  were  absent);  the  properties  of  these  robotic 
cilia  are  listed  in  table  A-5.  The  doublet  driving  mechanism  shown  in  figure  A- 12  was  used.  In 
figure  A- 15a.  the  lower  level  is  being  turned,  and  in  figure  A- 15b,  the  upper  lever  is  being 
turned.  In  both  cases,  the  cilium  is  bending  in  the  plane  normal  to  the  level  axis,  and  there  is  no 
torsion. 

Table  A-6  summarizes  the  in-plane  and  off-plane  curvature  tendencies  and  torque  required. 
In  summary,  the  sheathing  and  hardness  of  A53-63  and  A20-30  prevent  manifestation  of  the 
torsion;  the  radius  of  the  curling  is  lower  in  A20-30  because  the  curvature  k  is  higher  (radius  of 
curvature  =  1  Ik),  less  torque  is  needed  to  turn  the  lever  the  same  angle;  and  torsion  is 
manifested  only  in  the  case  when  torsional  rigidity  is  nearly  zero.  The  hardness  needs  to  be 
below  A20-30  for  the  torsion  to  be  manifested. 
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Figure  A- 1 5.  Cilia  Constructed  of  Continuous  Rubber 
(see  Table  AS)  (Hardness  A53-A63) 

(Broken  lines:  Lever  and  baseline  robotic  cilium  axes.) 


Table  A-5.  Physical  Properties  of  Two  Low-Viscosity  Two-Component  Urethane 
Casting  Elastomers  Used  for  the  Robotic  Cilia  Shown  in  Figure  A-14 
( Source :  H.  B.  Fuller  Company,  PO  Box  64683,  St.  Paul,  MN  55164) 


Property 

IJltralite  FH-3140 

Ultralite  FH-3138 

These  elastomers  are  designed  to 
make  flexible  molds  for  use  with 
thermosetting  resins;  they  pour 
and  cure  at  room  temperature. 

Used  to  make  tough,  yet  flexible 
molds;  not  moisture  sensitive 

Soft  and  flexible,  yet  tough 
compound  can  reproduce  finest 
detail  and  will  flex  to  release 
from  undercuts  w  ithout  tearing 

Shore  A  Hardness 

53-63 

20-30 

Specific  Gravity 

1.03 

0.992 
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Table  A-6.  Types  of  Distortion  Observed  in  a  Long  Robotic  Citium 
When  a  Torque  Is  Applied  at  the  Base 


Sample 

No. 

Hardness 

Length 

Observ  ation  of  Distortion 

W  hen  Moment  Is  Applied  at 
the  Base 

Torque  Required 
to  Shear  1°  per 

1  cm  of  Axial 
Length;  Scale:  0 
(low),  10  (high) 
(not  yet  measured) 

1 

Little  resistance  to  torque/deg 
due  to  the  presence  of  gaps 
between  segments 

Full 

Perfect  helical  spiral  of 
uniform  diameter  and  pitch 
(figure  A- 13) 

1-3 

2 

Little  resistance  to  torque/deg 
due  to  the  presence  of  gaps 
between  segments 

Half 

Spiral 

2 

2 

Sheathed  sample  2  resulting  in 
high  hardness 

Half 

Planar  bending  only;  no 
spiraling  (figure  A- 14) 

9-10 

3 

A20-30 

Half 

Planar  bending  only;  no 
spiraling 

4-5 

4 

A53-63 

Half 

Planar  bending  only;  no 
spiraling  (figure  A- 14) 

9-10 

A.6  PROGRESS  TOWARD  ROBOTIC  EMULATION  OF  NATURAL  SWIMMING 

Dynamic  systems  modeling  (see  figure  2  in  the  main  text)  has  shown  that  the  biological 
ciliurn  track  is  a  limit  cycle.  In  this  section,  the  effects  of  parameters  on  the  quality  of  agreement 
of  the  model  to  the  measurements  of  the  biological  ciliurn  are  considered,  along  with  the 
progress  made  in  the  robotic  emulation  of  natural  swimming. 


A.6.  /  Effects  of  Equation  Parameters  on  the  Dynamic  Systems  Modeling  Results 

The  exact  form  of  the  limit  cycle  in  equation  (1)  depends  on  the  values  of  the  main 
parameter  a  (see  Bandyopadhyay  et  al.,  2008).  The  effects  of  a  on  the  comparison  between  the 
Lienard  equation  model  (equation  (1))  and  measurements  is  shown  in  figure  A-16.  The  units 
shown  are  arbitrary:  the  modeling  results  are  in  arbitrary  units  of  volts;  and  the  units  (pm)  of  the 
biological  ciliurn  tracks  have  been  removed  to  scale  them  to  the  model  (the  v  results  are  not 
shown). 
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03 


Zvs  W  vs  V  for  9=0  021 


08 


Zdot  vs  Wdot  vs  Vdot  for  a=0  021 


Figure  A- 16.  Effects  of  the  Parameter  a  on  Olivo-Cerebellar  Modeling  of  Cilium  Motion 
(Model:  blue;  cilium  distal  point:  red.  Left  column:  position.  Right  column:  velocity. 
Values  of  a  are:  (a)  0.021,  (b)  0.022,  and  ( c)  0.030.) 
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A.6.2  Electronic  Analog  Dynamic  Systems  Actuation  of  Robotic  Cilium 

To  provide  theoretical  support  to  the  patent  applications  on  the  nonlinear  transduction 
described  in  Bandyopadhyay  (201 1  and  2012b),  this  section  describes  how  the  robotic  cilium  can 
be  actuated  using  the  dynamic  systems  model  in  equation  (1).  Figure  A-l  7  shows  the  apparatus 
for  future  use  in  documenting  the  perturbation  sensitivity  of  a  cilium  (disturbance  rejection)  to 
which  torque  is  applied  and  which  takes  a  limit  cycle  track  (see  figure  2).  (Also  see 
Bandyopadhyay  et  al.,  2008.) 


IO  Board 


Ethel  net 


IO  Host  j 


rtui 

D  VQ  Ctrl  | 


DAQ 


Position  C" 


Figure  A-l  7.  Apparatus  for  Driving  the  Robotic  Cilium  Using  the 
Electronic  Rendition  of  Dynamic  Systems  Equation  ( 1 ) 

The  schematic  on  the  right  in  figure  A-l  7  describes  the  relationship  of  the  components  in  the 
upper  left  and  lower  left  photographs.  The  lower  left  photograph  is  from  figure  1 1  in  the  main 
text  and  shows  the  robotic  cilium  and  the  four-motor  IOOM  drive.  The  upper  left  photograph 
shows  the  analog  circuit  (labeled  IO  for  inferior-olive  motion  control  neuron  (Llinas  et  al., 

2004))  for  generating  the  solution  to  equation  (1),  the  digital  controller,  and  the  digital  data 
acquisition  circuit  (DAQ:  digital  to  analog  and  analog  to  digital)  for  relating  the  analog 
component  to  the  digital  IOOMs.  The  wiring  is  shown  in  the  schematic  on  the  right.  (A  video 
show  ing  the  operation  of  the  robotic  cilium  using  the  analog  (IO)  device  is  available  from  the 
lead  author.) 
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A.7  SYNTHETIC  MODELING  OF  THE  TRANSFORMATION  FROM  SMALL  TO 

LARGE  VIA  OPTIMIZATION  MECHANISMS  OF  METACHRONISM 

AND  SELF-REGULATION 

#  'J 

A  meta  analysis  by  Bandyopadhyay  (201 1)  shows  that  the  propulsion  density  (kW/m)  of 
small  and  large  swimming  animals  (considering  their  red  muscles  only  for  cruising)  remains 
invariant  of  whether  they  are  large  or  small.  (In  fact,  this  relationship  extends  to  engineering 
underwater  swimming  vehicles,  irrespective  of  the  propulsion  mechanism.) 

Is  there  any  relationship  between  the  microscopic  cilium  propulsors  of  the  paramecium 
(Reynolds  number  Re  ~  0.1)  and  the  larger  propulsors  of  fish  (Re  ~  1  to  10  to  100's)  (which  have 
flexible  pectoral  flapping  fins)  and  even  larger  animals  such  as  penguins,  dolphins,  and  whales 
(which  have  larger,  stiffer  fins)  (Re  >  1000)? 

Information  on  how  the  oceans  changed  over  time  is  contained  to  some  extent  in  the 
boundary  conditions  and  optimization  rules  followed  by  swimming  animals.  To  aid  the 
hydrodynamic  study  of  such  synthetic  evolution,  the  relationship  of  cilium  actuator  to  eel  and 
fish  fins  was  notionally  explored  and  is  given  in  Bandyopadhyay  (2009). 

Consider  the  following  two  steps  showing  how  one  could  conceivably  go  from  small  (S)  to 
medium  (M)  and  from  small  to  large  (L),  as  shown  schematically  in  figures  A-l  8  through  A-21 : 

1.  Metachronism  as  the  Mechanism  of  Optimization  for  Transformation  from  S  to  M. 

Apply  metachronism  as  the  dominant  mechanism  to  optimize  from  S  to  M,  the  Reynolds  number 
range  being  smaller  than  that  in  S  to  L.  Assume  that  S  to  M  is  a  “bifurcation”  whereby  higher 
orders  of  freedom  are  introduced  while  retaining  the  intrinsic  dynamic  system  (equation  (1)). 

The  metachronic  coupling  of  the  downstream  cilia  may  now  be  added  to  generate  the  sinuous 
motion  of  the  dorsal  fin.  Metachronism  is  a  two-dimensional  mechanism  and  is  applied  to  the 
motion  of  a  two-dimensional  array  of  cilia,  shown  in  figure  A-l  9.  The  cilium  motion  is  three- 
dimensional,  and  the  array  is  wrapped  longitudinally  with  a  membrane,  whereby  spanwise  eel 
motion  can  be  produced.  Time  increases  from  (a)  to  (d);  the  bunching  of  the  cilia  denotes  a 
switch-over  from  the  power  stroke  to  the  return  stroke.  The  inset  in  figure  A- 19  is  a  top  view 
showing  the  D-tracks  of  the  cilia  and  how  they  alternate  in  the  transverse  plane  with  the  fin 
wave.  Qualitatively,  the  synthetic  trans-species  evolution  result  is  a  morphological  motion  of  the 
eel  dorsal  fin  from  the  cilium  motion.  This  hypothesis  requires  verification. 

2.  Self-Regulation  as  the  Mechanism  of  Optimization  for  Transformation  from  S  to  L. 

Assume  that  self-regulation  continues  as  Re  is  increased  from  S  to  L  and  that  equation  (1) 
applies.  The  re-mix  of  the  following  variables  is  considered  in  this  rule:  material  (G,  E)  and 
conformational  properties  (/),  Re  range,  lowering  of  the  natural  frequency  with  the  increase  in  Re 
and  forcing  at  this  frequency,  and  the  parameter  a  in  equation  (1). 

Consider  the  progression  from  S  to  L  directly.  The  cilium  has  nine  pairs  of  push-pull 
microtubules  spread  over  a  circular  rim  and  two  microtubules  at  the  center  for  actuation  (figure 
A-20a).  To  generate  L,  imagine  that  the  cilia  are  unwrapped,  as  shown  in  figures  A-20b  and 
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Figure  A-18.  Thrust  Actuators  from  S  to  M  to  L  Considered  in  Bandyopadhyay  (2009) 


Dorsal  fin 

(Membrane  over  cilia  array) 

/ 


(c)  (d) 


Figure  A-19.  Schematic  of  How  a  Metachronic  Wave  Produced  with  a  Densely  Packed 
and  Membrane-Wrapped  Two-Dimensional  Array  of  Cilia  Can  Emulate  an  Eel 

(Time  increases  from  a  to  d.) 
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(C)  (d) 

Figure  A-20.  Schematic  of  Progressive  Unwrapping  of 
the  Cilium  to  Produce  a  Fish ’s  Pectoral  Fin 


A-20c.  Assume  that  as  in  the  cilium  (Dute  and  Kung,  1978),  the  fin  in  L  also  has  a  virtual  origin, 
and  the  push-pull  microtubules  are  spread  out  from  their  virtual  origin  as  a  fan  and  are  shortened 
in  the  sequence  shown  in  figure  A-20d.  The  end  result  is  a  fish's  pectoral  fin  (see  Lauder  and 
Drucker,  2004).  Each  push-pull  pair  continues  to  be  welded  at  the  tip  (hemitrich),  whereby  a 
sliding  motion  at  the  base  would  bend  the  fin  bone  as  shown.  To  generate  larger  forces  at  higher 
Re.  the  fins  are  sheathed  as  in  the  eel  optimization  above.  Oscillatory  parameters  are  determined 
by  minimization  of  drag. 

The  fin  sequence  in  figure  A-21a  is  from  Drucker  and  Lauder  (2004),  with  a  virtual  origin 
added.  Figure  A-21d  shows  the  fin  undergoing  a  beat  cycle  similar  to  the  motion  of  a  fish's 
pectoral  fin,  as  in  the  video  of  Lauder  (2009). 


A-30 


9 

Virtual  Origin 


(a) 


Figure  A-21.  Sequence  Rationalizing  the  S  to  L  Progression 
(Arrows  indicate  the  direction  of fin  motion.  Oscillation  of  the  fin  assembly 
about  the  virtual  origin  qualitatively  reproduces  the  motion  of  the  sunfish 
pectoral  fin  as  in  Lauder  (2009).) 
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